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SUMMARY 

Stress  concentrations  near  grain  boundaries,  precipitates,  and  similar  micro-heterogeneities  nucleate 
instabilities  leading  to  macroscale  fracture.  As  it  is  not  practical  to  model  each  flaw  explicitly,  their  ensem¬ 
ble  effect  is  modeled  statistically.  Accounting  for  this  aleatory  uncertainty  requires  smaller  specimens  (e.g., 
small  finite  elements)  to  have  generally  higher  and  more  variable  strengths,  which  is  necessary  for  the 
initial  failure  probability  of  a  finite  domain  to  be  unaffected  by  its  discretization  into  elements.  Localiza¬ 
tion  itself,  which  might  be  attributed  to  constitutive  instability,  requires  realistic  numerical  perturbations 
to  predict  bifurcations  such  as  radial  cracking  in  axisymmetric  problems.  These  perturbations,  stemming 
from  microscale  heterogeneity,  are  incorporated  in  simulations  by  imposing  statistical  spatial  variability 
in  the  parameters  of  an  otherwise  conventional  (deterministic  and  scale-independent)  damage  model.  This 
approach  is  attractive  for  its  algorithmic  simplicity  and  straightforward  calibration  from  standard  strength 
tests.  In  addition,  it  results  in  virtually  no  loss  of  efficiency  or  robustness  relative  to  deterministic  models 
and  accommodates  general  three-dimensional  loading.  Despite  these  advantages,  some  significant  chal¬ 
lenges  remain  and  are  discussed.  However,  it  is  demonstrated  that  including  aleatory  uncertainty  with 
associated  scale  effects  significantly  improves  predictiveness  on  large-scale  computational  domains,  where 
it  is  impractical  to  resolve  each  crack  or  localization  zone.  Copyright  ©  2014  John  Wiley  &  Sons,  Ltd. 
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1.  INTRODUCTION 

This  paper  has  been  submitted  as  part  of  a  special  issue  in  honor  of  Professor  Ted  Belytschko, 
whose  contributions  are  pervasive  in  the  finite-element  community.  As  an  example,  the  production 
finite-element  code  that  we  use  in  this  paper  [1]  is  based  on  his  seminal  publications  on  hexa- 
hedral  elements  [2,  3].  Among  his  numerous  contributions  in  the  area  of  fracture  mechanics,  the 
extended  FEM  [4]  has  enabled  accurate  modeling  of  discontinuities  in  finite-element  meshes  (often 
introduced  by  cracks).  As  computational  resources  continue  to  evolve,  methods  such  as  extended 
FEM  might  ultimately  be  used  ubiquitously  to  model  failure.  Meanwhile,  damage  models,  such  as 
the  one  summarized  in  this  paper,  remain  useful  for  modeling  ballistic  impact  and  similar  high- 
energy  events.  We  have  selected  this  topic  for  this  special  issue  because  of  Professor  Belytschko’s 
continual  support  of  our  efforts  to  include  strength  statistics  and  scale  effects  in  engineering  damage 
simulations. 

This  paper  describes  a  practical  method  for  extending  a  computational  framework  to  accommo¬ 
date  aleatory  uncertainty  and  associated  size  effects  (e.g.,  spatial  variability  in  strength  and  other 
material  properties)  without  prohibitively  increasing  computational  overhead  or  model  development 
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Figure  1.  Tungsten-carbide  (WC)  sphere  impacting  a  silicon  carbide  (SiC-N)  ceramic  at  375  m/s.  All 
simulations  use  the  same  constitutive  model  on  the  same  three  mesh  resolutions,  but  those  in  the  bottom  row 
are  initialized  using  experimental  data  for  aleatory  uncertainty  and  associated  size  effects,  giving  relatively 
mesh-insensitive  radial  cracking  patterns  similar  to  those  in  the  (truly  representative)  experiment  [7]. 


time  relative  to  “conventional”  (deterministic  and  scale-independent)  damage  models.  As  illustrated 
in  Figure  1,  a  conventional  model’s  predictions  can  be  made  more  realistic  and  less  mesh  sensitive 
when  it  is  run  using  properties  that  include  statistical  variability  and  scale  effects  that  are  evident  in 
standard  laboratory  strength  testing.  Heterogeneity  is  applied  to  an  existing  conventional  smeared- 
damage  model  at  the  host-code  level  so  that  few  (if  any)  alterations  to  the  underlying  material 
model’s  source  code  are  needed.  Merits  and  drawbacks  of  the  conventional  damage  model  used 
for  this  study  [5,  6]  receive  only  slight  attention  in  this  paper,  whereas  the  physical  foundations, 
algorithms,  and  calibration  procedures  for  uncertainty  and  size  effect  are  covered  in  detail. 

Aleatory  uncertainty  refers  to  natural  spatial  variation  in  flaw  or  inclusion  morphology  that  causes 
local  regions  (at  a  scale  far  smaller  than  a  finite  element^)  to  be  weaker  or  stronger  than  surrounding 
material.  Because  it  is  impractical  to  precisely  measure  and  model  each  flaw  size  and  orientation, 
certain  material  properties  such  as  strength  are  distributed  statistically  in  a  manner  consistent  with 
macroscale  observations  for  a  variety  of  specimen  sizes.  Slight  spatial  variability,  which  is  relatively 
inconsequential  in  stable  deformations  [9, 10],  activates  localized  zones  of  damage  (e.g.,  shear  bands 
or  fractures)  as  the  intensity  of  loading  approaches  an  unstable  state. 

Conventional  plasticity  and  damage  models  have  long  been  known  to  admit  bifurcation  to  a 
localized  state  (when,  e.g.,  the  acoustic  tensor  has  a  zero  eigenvalue;  cf.  [12,  13]).  In  the  laboratory, 
however,  damage  begins  at  stress  levels  well  below  this  mathematical  stability  limit  [14,  15].  The 
situation  is  analogous  to  buckling:  a  column  will  buckle  at  a  load  well  below  the  theoretical  limit 
load  if  its  heterogeneities  (holes  and  pits)  are  large  enough  to  induce  non-infinitesimal  perturbations 
in  the  stress  field.  If  mesh  morphology,  material  properties,  and  loading  in  a  simulation  are  perfectly 
homogeneous,  then  conventional  deterministic  damage  models  (including  those  with  bifurcation 
and  decohesion  criteria)  are  mathematically  incapable  of  predicting  heterogeneous  fragmentation 
of  the  type  depicted  in  Figure  2.  Without  perturbations,  an  idealized  simulation  of  homogeneous 
loading  (i.e.,  run  at  “infinite  precision”)  would  reach  a  stability  limit  at  the  same  moment  for  all 
elements,  producing  unrealistic  failure  in  every  element  in  a  computational  domain,  which  is  an 
undesired  material  response  that  has  thwarted  numerous  attempts  to  simulate  dynamic  failure  of 
armor  ceramics  [16,  17].  Because  constitutive  instabilities  appear  to  be  extremely  sensitive  to  the 
nature  of  perturbations,  realistic  perturbation  fields  must  be  used  rather  than  relying  on  spurious 
mesh-dependent  and  platform-dependent  disturbances  such  as  numerical  round-off  (cf.  [18,  19]). 


^Simulations  presented  in  this  paper  were  run  using  a  produetion-quality  finite-element  eode  [1],  but  the  methods  apply 
equally  well  to  other  PDE  solvers.  The  same  techniques  have,  for  example,  been  confirmed  to  give  similar  results  when 
finite  elements  are  reinterpreted  as  particle  domains  in  a  material-point-method  code  [8]. 
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Figure  2.  Heterogeneous  dynamic  response  to  nominally  homogeneous  (and  initially  quasistatic)  unconfined 
compression  of  a  ceramic  (SiC-N).  The  simulation  uses  aleatory  uncertainty  in  strength  (i.e.,  macroscale 
statistical  variability  attributed  to  micro-heterogeneity)  to  capture  the  tendency  of  these  experiments  to  fail 
at  only  one  end,  typically  leaving  an  intact  fragment  at  the  other  end  [11]. 


(a)  Deterministic  model  (b)  Statistically  perturbed  and  scale-dependent  initialization  of  the  same  model 


Figure  3.  (a)  Spurious  mesh-texture-biased  radial  cracking  in  a  conventional  damage  model  compared  with 
(b)  realistic  randomly  oriented  cracking  (with  statistically  reproducible  spacing)  with  aleatory  uncertainty. 
For  a  SiC-N  cylinder  impacted  by  a  WC  sphere  at  375  m/s,  the  first  two  simulations  in  sub-figure  (b)  use  the 
same  mesh  as  sub-figure  (a) — each  with  a  different  random  seed,  whereas  the  third  uses  a  “paved”  mesh. 


A  damage  model  is  successful  if  it  is  able  to  reproduce  (without  mesh  sensitivity)  observed  scale- 
dependent  and  rate-dependent  statistical  bifurcation  into  damaged  and  intact  zones  with  realistic 
expended  fragmentation  energy. 

In  a  systematic  investigation  of  correlations  between  standard  material  properties  (such  as 
hardness,  grain  size,  porosity,  etc.)  and  ballistic  performance  of  a  variety  of  SiC-N  ceramics,  Ray 
et  al.  [20]  found  that  “a  wide  variety  of  SiC-based  materials  can  give  good  ballistic  results,  con¬ 
tradicting  some  of  the  theories  about  what  is  important  to  improve  ceramics  for  armor.”  Their 
investigation  revealed  a  possible  correlation  between  the  Weibull  modulus  and  ballistic  performance, 
but  they  were  quick  to  note  that  this  result  might  have  been  coincidental  because  there  was  also 
no  correlation  to  the  weak  tail  of  the  strength  distribution.  If  these  conclusions  are  independently 
validated  through  further  testing,  then  variability  itself  (not  the  value  of  the  lowest  strength)  might 
be  the  key  factor  in  ballistic  performance,  suggesting  that  the  salient  numerical  requirement  is  to 
realistically  seed  an  instability,  perhaps  of  the  nature  discussed  by  Grinfeld  [21]. 

For  axisymmetric  loading,  such  as  indentation  of  a  macroscopically  homogeneous  cylinder, 
symmetry-breaking  radial  cracks  originate  from  micromorphological  instability  [21].  Just  as  a 
mathematical  perturbation-free  analysis  of  column  buckling  produces  unrealistic  simultaneous 
failure  at  every  point  in  the  column,  a  perturbation-free  mathematical  analysis  of  axisymmetric 
loading  produces  simultaneous  failure  of  all  points  equidistant  from  the  symmetry  axis  or  excessive 
cracking  in  a  deterministic  simulation  such  as  Figure  3(a).  Round-off  noise  represents  a  numeri¬ 
cally  infinitesimal  perturbation  many  orders  of  magnitude  smaller  than  the  perturbations  evident 
in  experimental  data.  As  illustrated  in  Figure  3,  measured  statistics  in  strength  should  be  used  to 
establish  heterogeneous  material  property  fields  of  statistically  homogeneous  media.  When  simula¬ 
tions  rely  exclusively  on  round-off  error  to  stimulate  instabilities,  the  results  might  break  symmetry 
without  explicit  control  of  perturbations,  but  such  predictions  have  no  physical  basis,  as  they  are  typ¬ 
ically  caused  by  the  following:  (1)  mesh-dependent  or  platform-dependent  differences  in  the  fields 
or  (2)  order  of  operation  errors  in  the  code  (where,  for  example  the  stress  field  at  certain  locations 
is  updated  at  the  wrong  time).  In  Figure  3(b),  realistic  random  crack  orientations  exhibit  reduced 
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texture  bias  on  the  radial-trisection  meshes.  Timmel  [22]  notes  that  other  mesh  types  might  pro¬ 
mote  realistic  failure  patterns,  but  algorithms  that  minimize  this  mesh-texture  bias  are  preferable. 
The  “paved”  mesh  in  the  rightmost  image  in  Figure  3(b),  for  example,  encourages  slight  curva¬ 
ture  in  the  crack  paths  without  any  other  pronounced  differences  in  comparison  with  the  radial 
trisection  results. 

In  addition  to  needing  realistic  statistical  strength  perturbations  in  simulations,  an  equally 
important  feature  is  the  well-known  concomitant  dependence  of  strength  on  specimen  size  [23].  As 
material  failure  tends  to  originate  at  the  weakest  point  in  a  material,  large  samples  will  be  weaker 
(on  average)  because  they  are  more  likely  to  contain  a  critically  large  or  critically  oriented  flaw. 
Because  the  location  of  the  weakest  points  in  a  material  may  not  coincide  with  the  points  of  high¬ 
est  stress,  localized  material  failure  does  not  generally  originate  at  the  highest  stress  concentrations, 
nor  does  it  necessarily  precipitate  cascading  formation  of  macrocracks.  With  aleatory  uncertainty 
in  simulations,  non-cascading  local  damage  (i.e.,  non-weakest-link  behavior)  is  evident  in  the 
“speckled”  failure  points  in  Figure  3(b)  (which  are  very  similar  to  those  observed  in  computed 
tomography  data;  cf.  [24],  which  likewise  do  not  all  coalesce  into  macrocracks). 

Using  statistics  in  constitutive  modeling  dates  back  to  the  19th  century,  where  work  was  driven  by 
geomaterials  because  of  their  heterogeneity  even  at  the  macroscale.  Statistics  have  been  also  used 
extensively  to  predict  concrete  fracture  [25]  and  munitions  fragmentation  or  “hot-spot”  formation 
in  energetic  media  [26-28].  Convergence  of  the  statistical  distribution  of  fragments  generated  by  a 
uniformly  expanding  ring  was  demonstrated  by  including  Weibull  variability  of  the  yield  strength 
(without  a  size  effect)  [29].  Weibull  statistics  and  size  effects  applied  to  cohesive  zone  models 
improve  their  numerical  convergence  and  predictive  properties  [30],  without  necessarily  implying 
that  the  overall  strength  is  itself  Weibull  distributed. 

In  Section  2,  we  summarize  a  recent  study  [31]  showing  that  strength  distributions  strongly 
depend  on  the  tensorial  character  of  the  loading  mode,  which  leads  to  the  mathematically  rigorous 
concept  in  Section  3  of  “quantile”  (isoprobability)  surfaces  that  generalize  the  limit  functions 
of  deterministic  damage  theories.  Stress  dependence  of  failure  is  characterized  in  Section  4  and 
Section  5  from  examination  of  experimental  data,  which  suggest  that  both  scale  effects  and  aleatory 
uncertainty  decrease  with  increasing  pressure.  These  observations  motivate  using  a  Weibull  strength 
distribution  only  for  spherical  (isotropic)  tension,  whereas  non-Weibull  distributions  (as  well  as 
decreased  statistical  variability  with  pressure)  are  naturally  predicted  for  states  having  a  nonzero 
stress  deviator.  Data  delocalization  methods  accounting  for  non-uniform  stress  fields  in  labora¬ 
tory  testing  are  discussed  in  Section  6,  along  with  recent  work  on  so-called  data  reloealization 
needed  in  simulations  to  compensate  for  the  inability  of  low-order  shape  functions  to  describe  stress 
variation  on  an  element  [32].  The  constitutive  damage  model  used  in  our  case  studies  is  briefly  dis¬ 
cussed  in  Section  7  (though  our  approach  applies  to  any  damage  model).  Finally,  although  failure 
initiation  receives  the  most  attention  in  this  paper.  Section  8  emphasizes  the  need  for  including 
scale-dependent  damage  rate  effects  during  the  failure  progression  phase  of  the  simulation.  Failure 
initiation  (e.g.,  the  first  growth  of  a  microcrack  leading  to  initiation  of  stiffness  degradation)  is 
herein  treated  as  a  weakest  link  phenomenon,  but  there  is  no  assumption  that  initiation  of  failure 
will  necessarily  progress  to  a  complete  loss  of  stiffness  or  strength.  Thus,  overall  structural  failure 
is  not  generally  a  weakest  link  event. 

2.  EXTENSION  OF  WEIBULL  THEORY  FOR  GENERALIZED  STRESS  STATES 

Weakest  link  theory  (which  we  do  not  assume  applies  to  all  loading  modes)  is  governed  by  Weibull 
strength  distributions,  leading  to  Weibull  scale  effects  in  which  large  specimens  (and  hence  large 
finite  elements)  are  weaker,  on  average,  than  small  specimens  (or  small  finite  elements).  For  general 
loading  trajectories,  our  phenomenological  extension  of  Weibull  theory  will  be  shown  to  exhibit 
realistic  non-Weibull  statistical  strength  distributions  and  non-Weibull  scale  effects. 

2.7.  Failure  theories  and  implications  for  specimen  size  dependence 

If  it  were  possible  to  create  a  sample  containing  exactly  one  small  flaw,  then  onset  of  failure  under  a 
given  tensor  stress  state  T  (depicted  in  Figure  4  as  a  3D  Mohr’s  circle)  relates  to  the  likelihood  that 
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Figure  4.  Uncertainty  in  a  single  crack’s  orientation  (quantified  by  its  unit  normal,  which  is  a  point  on  the 
unit  sphere)  allows  the  stress  state  (Mohr’s  circle)  to  fall  outside  a  failure  line,  depicted  here  as  a  straight 
line  corresponding  to  a  criterion  that  the  crack  will  fail  if  the  resolved  shear  stress  (minus  friction)  reaches  a 
critical  value.  Statistically  variable  crack  sizes  and  failure  modes  perturb  the  failure  line  itself  (cf.  [31]). 

the  flaw  is  critically  oriented  or  critically  large.  Figure  4,  for  example,  shows  that  even  a  simplistic 
Mohr-Coulomb  failure  criterion  has  uncertainty  of  failure  if  the  crack  orientation  is  random  [31]. 

Similar  ohservations  have  been  made  using  more  sophisticated  single-crack  failure  theories  [33], 
and  the  same  concepts  have  been  applied  to  account  for  variable  slip-plane  orientations  in  plasticity 
[34].  Some  flaw  orientations  are  safe  from  failure,  whereas  others  (dark  regions  on  the  spheres  in 
Figure  4)  are  not.  Brannon  and  Gowen  [31]  limit  their  scope  to  failure  initiation  in  a  statistically 
isotropic  medium  (for  which  flaw  orientation  is  initially  uniformly  random  so  that  the  initial  failure 
probability  is  the  ratio  of  solid  angle  of  the  zone  of  critical  orientations  to  the  area  of  the  sphere). 
Using  analytical  methods,  they  derive  failure  statistics  (probabilities  that  samples  will  fail  under  spe¬ 
cific  loading  conditions)  from  uncertainties  in  flaw  orientations,  and  they  illustrate  the  corresponding 
concept  of  strength  quantile  surfaces  that  generalize  the  notion  of  a  limit  failure  surface  used  in 
classical  damage  theories. 

Dienes  et  al.  [28]  and  Nemeth  [35]  accommodate  intrinsic  and  deformation-induced  anisotropy  by 
allowing  flaw  size  distributions  to  be  correlated  to  flaw  orientation  either  from  the  outset  or  through 
oriented  crack  growth.  Macon  et  al.  [36]  further  point  out  that  induced  anisotropy  (leading  to  what 
Nemeth  [35]  refers  to  as  “equatorial”  flaw  distributions)  may  represent  a  way  to  model  dilatation  in 
triaxial  compression  in  a  way  that  avoids  anomalous  predictions  of  hydrostatic  strengthening  under 
dilatation  seen  in  some  isotropic  hardening  models  (e.g.,  [37]).  Rather  than  further  expanding  on  the 
physical  motivations  for  these  so-called  unit-sphere  representations  of  failure  statistics,  we  focus  on 
actual  implementations  of  these  theories  in  engineering  simulations. 

2.2.  Weibull  theory  {and  its  extension)  for  generalized  stress  states 

Consider  a  stress  state  T  =  aN,  where  a  is  the  magnitude  of  stress  and  N  is  a  unit  tensor  called 
the  stress  trajectory  or  loading  mode.  For  example,  uniaxial  stress  has  N  =  DIAG[1, 0, 0],  whereas 
spherical  tension  has  N  =  DIAG[1, 1,  l]/\/3  and  biaxial  tension  has  N  =  DIAG[1, 1, 0]/a/2-  For 
simplicity,  our  discussions  might  presume  a  straight  path  through  stress  space  (where  N  is  constant), 
but  the  concept  of  strength  quantile  surfaces  in  six-dimensional  (6D)  stress  space  [31]  admits  more 
tortuous  tensorial  paths  through  stress  space  in  which  the  scalar  a  (called  strength)  is  merely  the 
value  of  the  path  parameter  at  failure.  Just  as  deterministic  plasticity  yield  surfaces  generally  evolve 
with  plastic  flow,  strength  quantile  surfaces  are  also  expected  to  change  in  a  history-dependent 
manner  so  that  probabilistic  strength,  a,  generally  depends  implicitly  on  evolving  internal  state  vari¬ 
ables  such  as  crack  density  in  addition  to  the  loading  mode.  Let  a  be  the  median  strength  observed 
when  testing  a  sample  of  volume  V  under  a  given  loading  mode  N.  To  account  for  microscale  het¬ 
erogeneities  (e.g.,  stress  concentrations  at  flaws)  and  to  accommodate  associated  scale  effects,  a 
theory  is  required  for  predicting  the  probability,  V),  that  a  sample  of  a  volume  V  (gen¬ 
erally  distinct  from  V)  will  begin  failing  under  an  applied  stress  T.  The  V)  function  is 

essentially  a  cumulative  distribution  function  (CDF)  in  6D  stress  space,  so  it  is  a  Heaviside  step 
function  in  a  deterministic  theory  (i.e.,  its  gradient  is  singular  at  the  limit  surface).  Otherwise,  for 
non-deterministic  theories,  the  directed  derivative  along  a  stress  trajectory  gives  the  probability 
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Figure  5.  CDFs  for  normalized  failure  stress  E  for  various  normalized  sample  volumes,  a,  and  various 

values  of  the  Weibull  modulus  m. 


density  function  (PDF)  for  trajectory-dependent  strength  as  p{a,  V)  =  (9P^“'/9T)  :  N,  which  is  the 
key  observation  needed  to  construct  the  physically  meaningful  tensor-hased  function,  V), 

whose  level  sets  are  dubhed  as  “aleatory  quantile  surfaces”. 

For  a  fixed  stress  trajectory  N,  define  non-dimensionalized  stress  and  size  ratios: 

a  V 

T,  —  —  and  a  —  =.  (1) 

a  V 

Weibull  theory  (which  we  soon  argue  does  not  apply  to  general  loading  modes  but  is  nevertheless 
useful  for  illustrating  concepts  in  a  familiar  context)  predicts  that  a  sample  of  non-dimensionalized 
volume  a  will  begin  failing  under  non-dimensionalized  stress  S  with  probability 

pfail  ^  j  _  psafe  ^  j  ^  (2) 

which  is  the  initial  CDF  for  strength,  in  which  m  is  a  constant  called  the  Weibull  modulus  (or 
shape  parameter).  The  complementary  CDF  psaf®(T,  V)  represents  the  probability  that  a  sample  of 
volume  V  will  not  have  begun  failing  under  an  applied  stress  T.  Plots  of  Equation  (2)  are  shown  in 
Figure  5  for  various  relative  sample  volumes  and  Weibull  moduli,  demonstrating  that  small  values 
of  m  correspond  to  large  variability  in  strength.  A  purely  deterministic  model  corresponds  to  a  unit 
step  function  in  the  limit  m  ^  oo.  Because  the  plots  for  small  specimens  (small  a)  fall  to  the  right 
of  those  for  larger  specimens,  small  specimens  have,  on  average,  higher  strengths.  Equivalently, 
because  curves  for  small  specimens  fall  below  those  of  larger  specimens,  a  small  specimen  is  less 
likely  to  fail  at  any  given  applied  stress.  Eurthermore,  the  plots  show  that  small  specimens  have 
greater  spread  in  their  strengths. 


2.3.  Mode,  median,  and  standard  deviation  of  the  Weibull  distribution 

Eor  the  Weibull  distribution  introduced  in  Equation  (2),  the  median  is  the  50%  quantile, 

representing  the  stress  level  at  which  half  of  the  specimens  in  a  large  population  will  have  initiated 
failure.  Hence,  S™®*™  is  the  value  of  S  at  which  =  1/2,  and  referring  to  Equation  (2),  it 
follows  that  the  median  strength  depends  on  specimen  size  according  to 

^median  _  1/m 

The  probability  density  function  (PDE),  />(S,  a)  is  the  derivative  of  the  CDE  with  respect  to  S: 

opfail 

p{i:,a)  ^  —  =  (In 2)  (2-“^“)  am  S'””!  (4) 

The  so-called  Weibull  scale  parameter  the  mode  (at  the  peak  of  the  PDF),  and  the 

mean,  S“®^"(a)  =  S  p{Yj,a)  dS  are  found  from 

j^scale  _  j^mode  f  _  \  \h  Smean  + 

- ^  =  (ln2)i^,  - -  ,  and  - ^ ^  (5) 

^median  '  ^  ^median  y  777  In  2  y  ^ median  [In  2]^'^ 

where  F  is  the  gamma  function.  Variability  (spread)  in  the  data  may  be  quantified  by  the  standard 
deviation,  ,  illustrated  in  Figure  6(a)  and  defined  by 
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Figure  6.  The  graphs  in  (a)  imply  that,  in  addition  to  higher  strength,  greater  variability  must  apply  at  smaller 
elements  within  a  computational  domain,  such  as  those  near  the  centerline  in  (b),  giving  an  unavoidable 
impression  of  heterogeneity  in  what  is  nevertheless  a  statistically  homogeneous  strength  distribution  so  that  a 
physical  domain’s  failure  initiation  probability  under  homogeneous  loading  is  unaffected  by  the  meshing  of 

that  domain. 


^median 


Jt 


s/6, 


for  large  m 


m 


(6) 


The  plot  of  the  standard  deviation  in  Figure  6(a)  illustrates  our  earlier  assertion  that  one 
should  not  rely  on  round-off  error  to  initiate  mathematical  instabilities.  Even  for  the  largest  Weibull 
modulus  in  Figure  6(a),  the  standard  deviation  in  strength  can  be  on  the  order  of  15%.  In  fact,  for  the 
left-hand  side  of  Equation  (6)  to  equal  a  typical  round-off  error  of  10“^®,  the  Weibull  modulus  would 
need  to  be  more  than  10  orders  of  magnitude  larger  than  reported  for  any  polycrystalline  material 
to  our  knowledge.  Hence,  perturbations  that  precipitate  instabilities  or  bifurcations  associated  with 
localization  in  a  numerical  code  must  be  realistically  seeded  in  order  to  obtain  realistic  results.  A 
deterministic  model  unrealistically  has  no  variability  (zero  standard  deviation)  in  strength. 

Hild  et  al.  [38]  identify  a  characteristic  sample  (element)  size  marking  a  statistical-deterministic 
transition  such  that  larger  elements  are  treated  deterministically  while  variability  is  applied  for 
smaller  elements.  That  approach  admits  nonphysical  jump  discontinuities  in  the  standard  deviation, 
whereas  Figure  6(a)  shows  that  seeding  all  elements  with  a  Weibull  distribution  results  in  a  smooth 
transition  toward  a  deterministic  response  (zero  standard  deviation)  as  the  element  size  increases. 
Both  higher  strength  and  greater  variability  are  evident  in  the  time-zero  initialization  of  a  statisti¬ 
cally  homogeneous  material  in  the  computational  domain  of  Figure  6(b),  whieh  was  generated  using 
concepts  developed  in  the  next  section  to  accommodate  smooth,  trajectory-dependent,  and  generally 
non- Weibull  seeding  of  strength  for  3D  stress  states. 

Failure  probability  functions  are  valuable  in  engineering  if  one  seeks  a  probability  of  failure 
initiation  for  the  purpose  of,  say,  component  reliability  and  safety  assessment  [39,  40].  Theories 
are  also  available  for  expected  values  of  fragment  sizes  in  impact  and  shock  loading  after  failure 
initiation  [41,  42].  However,  if  one  is  interested  in  detailed  analysis  of  the  dynamic  effects  of  the 
failure  process,  then  realizations  of  strength,  consistent  with  the  function,  are  required 

to  generate  realistic  crack  spacings  and  fragment  sizes  [43].  Simulating  realistic  damage  and  frag¬ 
mentation  events  is  especially  important  for  assessing  degradation  of  strength  that  could  lead  to 
particular  vulnerabilities  to  subsequent  impact  events.  The  next  section  of  this  paper  recognizes 
and  accounts  for  the  existence  of  a  y  s^  function  without  requiring  that  it  ever  be  explicitly 

constructed  in  simulations.  Accordingly,  we  will  now  focus  on  the  goal  of  converting  failure 
probabilities  into  strength  realizations  on  a  computational  domain. 
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3.  “FUZZY  FAILURE”  THEORY  AND  STRENGTH  REALIZATIONS 

Classical  deterministic  plasticity  models  presume  existence  of  a  scalar-valued  yield  function, 
/(T,  U),  such  that  a  stress  state  T  applied  to  a  specimen  of  volume  V  is  elastic  if  /(T,  U)  <  0.  The 
yield  surface  is  the  set  of  T  for  which  / (T,  V)  —  0.  Scale-insensitive  models  omit  dependence  on 
V .  The  yield  function  also  implicitly  depends  on  internal  variables,  such  as  dislocation  density  and 
crack  texture  parameters,  which  can  evolve  to  allow  the  yield  surface  to  expand  (harden)  or  contract 
(soften)  in  response  to  plastic  loading.  The  yield  surface  hounds  the  domain  of  elastically  attainable 
states;  engineering  damage  theories  also  use  a  so-called  limit  function  F(T,  V)  that  defines  the 
boundary  of  stress  states  attainable  by  any  quasistatic  means  [44].  A  quasi-brittle  model  employs 
distinct  yield  and  limit  functions  (/  ^  F),  whereas  a  purely  brittle  model  uses  coincident  func¬ 
tions  (/  —  F).  Damage  corresponds  to  collapse  of  the  limit  surface  and  concomitant  degradation 
of  elastic  stiffness. 

These  conventional  sharp-threshold  plasticity  models  have  required  considerable  time  and 
resources  to  develop.  Many  (such  as  [45]  and  [46])  are  also  equipped  with  large  material  property 
databases,  compiled  over  many  years.  By  presuming  that  existing  discrete  threshold  functions  are 
actually  fitted  only  to  median  experimental  data,  so-called  “fuzzy  failure  theory”  incorporates  spa¬ 
tial  variability  and  scale  effects  in  these  models.  Moreover,  because  of  the  size  effect,  these  fits  apply 
only  to  the  specimen  size  used  in  the  calibration  experiments,  where  specimen  size  is  defined  using 
“data-delocalization”  methods  as  explained  in  Section  6. 

5.7.  Aleatory  uncertainty  theory 

Conceptually,  spatial  variability  is  incorporated  into  these  conventional  plasticity  models  by 
replacing  a  limit  function  F(T,  V)  by  a  scale-dependent  probability  field  ranges 

continuously  from  a  value  of  1  at  stress  states  that  are  safe  from  failure  to  a  value  of  0  at  stress 
states  that  are  certain  to  induce  failure.  The  conventional  limit  function  could,  for  example,  be  iden¬ 
tified  with  the  50%  quantile  isosurface:  F(T,  V)  —  P'*af‘^(T,  U)  —  1/2.  Of  course,  similar  revision 
of  the  yield  function  applies  to  quasi-brittle  models.  Strength  is  statistically  variable  because  of 
aleatory  uncertainty  in  the  micromorphology  (such  as  crack  densities  and  textures).  Recalling  that 
F(T,  V)  depends  implicitly  on  internal  variables  that  directly  or  indirectly  account  for  micromor¬ 
phology,  it  follows  that  a  natural  way  to  indirectly  incorporate  the  F'*^^®(T,  V)  in  a  simulation  is 
to  retain  an  existing  code  infrastructure  for  a  scale-insensitive  F(T),  but  to  initialize  each  finite 
element  with  different  F(T)  functions  by  statistically  perturbing  the  parameters  defining  F(T)  in 
a  way  that  recovers  observed  strength  probabilities  captured  in  F'*^®(T,  V),  where  V  is  taken  as 
the  finite-element  volume.  As  this  initialization  phase  of  the  simulation  is  the  only  point  at  which 
psafe^^rp  needed,  there  is  no  algorithmic  burden  to  otherwise  incorporate  this  function  into 

run-time  calculations,  making  the  implementation  of  aleatory  uncertainty  quite  manageable. 

The  F'*^^®(T,  V)  function  allows  the  effects  of  micromorphology  to  be  treated  in  ensemble  through 
empirical  strength  CDFs.  With  this  engineering  approach,  the  median  (which  has  been  presumably 
already  fit  to  laboratory  data)  is  retained  but  scaled  appropriately  with  specimen  (finite-element) 
size.  Statistically  perturbed  realizations  of  strength  surfaces  in  stress  space  are  centered  about 
the  scale-adjusted  median  (i.e.,  50%  quantile)  limit  surface  according  to  the  statistics  evident  in 
laboratory  data  (and  consistent  with  limit  surface  admissibility  requirements  such  as  convexity). 

“Fuzziness”  in  the  limit  surfaces  of  conventional  plasticity  models  has  been  formalized  by 
Brannon  and  Gowen  [31],  who  introduce  the  concept  of  quantile  (isoprobability)  surfaces  in  stress 
space,  replacing  the  single  failure  surface  used  in  deterministic  models  by  the  family  of  surfaces 
in  Figure  7(a).  The  red  in  Figure  7(a)  depicts  the  median  (50%  quantile)  isoprobability  surface 
viewed  down  the  hydrostat,  with  the  5%  and  95%  ‘quantile  surfaces  depicted  as  blue.  The  images 
in  Figure  7(b)  are  similar  except  that  quantile  surface  colors  are  replaced  with  statistical  dot  densi¬ 
ties  centered  about  the  median  quantile  surface.  The  “weak  tail”  of  a  strength  distribution  refers  to 
the  smaller  quantile  surfaces,  which  (being  closer  to  the  origin  in  stress  space)  represent  instances 
of  low-stress  failure.  For  the  idealized  micromorphology  used  to  generate  Figure  7(a),  the  median 
50%  quantile  surface  (colored  red)  is  a  rounded  triangle  (modeled  using  the  Willam-Wamke  third- 
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Mohr-Coulomb  Willam-Warnke 


(a)  Quantile  surfaces 


(b)  Fuzzy  failure  surfaces  (octahedral  view) 


Figure  7.  Accounting  for  variability  in  flaw  size  affects  the  location  of  the  failure  threshold  itself,  producing 
quantile  surfaces  in  stress  space  depicted  here  looking  down  the  hydrostat  and  visualized  according  to  initial 
failure  probability  density  by  (a)  colors  showing  theoretical  shapes  of  each  quantile  surface  [31]  or  (b)  dot 
densities  that  were  generated  using  self-similar  quantile  shapes  (more  like  the  current  implementation). 


invariant  function  common  in  geomechanics;  cf.  [6]),  whereas  the  weak  5%  quantile  is  the  blue 
distorted  hexagon  that  is  easily  modeled  with  Mohr-Coulomb  third-invariant  dependence.  The 
quantile  functions  in  Figure  7(b)  did  not  include  this  recent  finding  of  Lode-angle  sensitivity. 

In  our  fuzzy  failure  theory,  natural  subscale  heterogeneities  are  accommodated  in  a  macroscale 
simulation  by  using  different  realizations  of  strength  in  each  finite  element,  such  that  the  median 
strength  of  a  finite  volume  is  preserved,  with  a  key  goal  being  that  failure-initiation  statistics  ought 
to  be  unaffected  by  mesh  refinement.  In  numerical  implementations,  each  finite  element  is  regarded 
as  a  finite-sized  specimen  requiring  assignment  of  strength  consistent  with  scale  effects  observed  in 
the  laboratory.  For  illustration  purposes,  suppose  that  the  stress  at  failure  a  for  a  specimen  size  V 
is  known  to  be  Weibull  distributed  with  median  a  for  a  given  stress  trajectory  N  (keep  in  mind  that 
non-Weibull  strengths  are  still  possible  for  other  trajectories).  To  assign  strength  properties  to  a  finite 
element,  one  must  generate  a  uniformly  distributed  random  number  R  on  the  interval  0  ^  i?  <  1 . 
The  realization  of  failure  stress  can  then  be  found  by  replacing  with  R  in  Equation  (2)  and 
solving  for  S.  The  result  is 


InR 


1/m 


a  ln(l/2) 


i.e.. 


a 


VlnR 

Fln(l/2) 


1/m 


(7) 


For  a  different  sample  volume  {V  ^  V),  this  equation  predicts  a  different  strength.  In  our  case,  the 
strength  a  assigned  to  a  finite  element  depends  on  the  volume  V  of  the  element  relative  to  the  sample 
volume  used  in  laboratory  testing.  Equal-sized  regions  always  have  the  same  probability  of  failure 
initiation  under  homogeneous  loading  regardless  of  whether  or  not  those  regions  are  subdivided 
into  many  or  few  elements.  As  a  special  case,  Equation  (7)  ensures  that  a  collection  of  small  ele¬ 
ments  having  total  volume  V  will  have  the  same  failure  initiation  statistics  (under  uniform  loading) 
as  observed  in  the  original  laboratory  sample  of  volume  V  (non-uniform  loading  is  addressed  in 
Section  6).  Equation  (7)  therefore  contains  an  implicit  length  scale  that  ensures  mesh-insensitive 
predictions  for  the  onset  of  failure  under  homogeneous  loading.  This  is  verified  numerically  for 
a  spherical  stress  trajectory,  N  =  DIAG[1, 1,  l]/v/3,  in  Section  5  (Eigure  14).  As  discussed  in 
Section  7,  mesh  insensitivity  for  the  subsequent  progression  of  failure  is  an  active  area  of  research  in 
which  nonlocal  effects  and  energy  release  considerations  are  important  for  mathematical,  numerical, 
and  physical  reasons. 

Plots  of  the  material  strength  associated  with  each  realization  parameter  R  in  Equation  (7),  and 
strength  distribution  curves  associated  with  multiple  realizations,  are  shown  in  Eigure  8. 


3.2.  Algorithm  for  implementing  aleatory  uncertainty 

A  computational  procedure  to  support  aleatory  uncertainty  (fuzzy  thresholds)  in  material  properties 
can  be  designed  such  that  an  existing  code  infrastructure  is  impacted  minimally  and  so  that  existing 
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Figure  8.  Strength  realization  functions  for  various  Weibull  moduli,  given  by  Equation  (7),  and  probability 
density  functions  for  failure  stress,  given  by  Equation  (4). 

conventional  deterministic  material  models  require  no  internal  revisions  of  their  subroutines.  Our 
approach  is  outlined  below 

(1)  Modify  the  user  input  parser  so  that  material  parameters  that  were  formerly  deterministic  and 
specihed  (for  example)  by  a  syntax 

PROP  =  value 

are  treated  statistically  whenever  the  user  employs  a  syntax  (for  example)  of  the  form 

PROP  =  Weibull  (median,  Vref,  Wmod) 

where  median  is  the  median  measured  in  the  laboratory  using  a  sample  of  volume  Vref  and 
Wmod  is  the  Weibull  modulus.  Because  median  usually  equals  value  from  the  deterministic 
model,  this  approach  to  aleatory  uncertainty  requires  two  new  user  inputs  (Vref  and  Wmod), 
both  of  which  will  be  available  if  the  original  strength  measurements  were  accompanied  by 
specimen  size  and  statistically  significant  standard  deviations  (as  should  be  the  case  if  testing 
followed  ASTM  and  similar  standards). 

(2)  For  each  statistical  user  input,  allocate  storage  for  one  field  variable,  giving  it  a  name  such  as 
PROPSR  (where  “PROP”  is  the  property  name  and  “SR”  stands  for  statistical  realization). 

(3)  Within  an  initialization  loop  over  each  finite  element,  compute  a  uniformly  distributed  random 
number  on  [0, 1),  and  use  the  finite-element  volume  V  in  Equation  (7)  (or  any  other  distri¬ 
bution  of  interest)  to  assign  a  realization  value  “a”  for  the  statistical  property.  Save  the  result 
in  the  PROPSR  field  variable  array. 

(4)  During  each  time  step,  before  calling  the  constitutive  model,  construct  a  copy  of  the  material 
property  array  in  which  PROP  is  replaced  by  PROPSR,  and  send  this  revised  property  array  to 
the  material  model. 

Virtually  all  of  the  CPU  overhead  resides  in  step  3,  which  is  a  once-only  initialization  setup  task 
typically  requiring  only  a  tiny  fraction  of  overall  CPU  for  the  calculation.  Because  step  4  requires 
nothing  more  than  a  small  amount  of  additional  memory  storage  and  an  array  copy,  this  approach 
to  aleatory  uncertainty  entails  almost  no  change  in  overall  CPU  costs  relative  to  the  deterministic 
model.  Because  the  constitutive  model  is  treated  as  a  “black  box”,  the  necessary  code  architecture 
may  be  written  in  a  way  that  supports  all  constitutive  models  in  the  code  and  not  just  the  damage 
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model  of  interest.  Users  must  be  cautioned,  however,  that  carelessly  or  arbitrarily  perturbing  param¬ 
eters  might  violate  the  model’s  admissibility  conditions  (such  as  the  convexity  requirements).  As 
discussed  in  the  next  section,  achieving  physically  reasonable  distributions  for  material  parameters 
might  be  facilitated  by  introducing  a  change  of  variables. 


4.  DIRECTIONALLY- WEIBULL  THREE-INVARIANT  THEORY 

The  aleatory  quantile  function,  probability  that  a  sample  of  size  V  is  safe  from 

failure  at  an  applied  stress  T.  This  function  is  used  indirectly  in  implementations  of  aleatory  uncer¬ 
tainty.  Specifically,  the  material  parameters  implicitly  contained  in  the  limit  function  F(T,  E)  are 
statistically  perturbed,  thus  giving  each  finite  element  its  own  perturbed  limit  surface  distinct  from 
other  elements.  If  F(T,  V)  depends  on  only  a  single  material  parameter  (as  for  von  Mises  or  Tresca 
theory),  then  model  calibration  requires  measured  strength  statistics  along  only  a  single  loading  tra¬ 
jectory  N  (such  as  simple  shear).  In  contrast  to  calibration,  validation  requires  additional  data  along 
other  trajectories  (such  as  uniaxial  compression,  biaxial  tension,  etc.).  Such  efforts  might  imply  the 
need  for  a  ffiree-parameter  limit  function,  for  which  calibration  requires  strength  data  along  three 
non-crossing  trajectories,  as  later  discussed  in  Section  5.2. 

Measured  strengths  are  interpreted  in  the  context  of  aleatory  quantile  (i.e.,  constant  failure  ini¬ 
tiation  probability)  surfaces  embedded  in  full  6D  tensorial  stress  space.  These  multi-dimensional 
surfaces  are  typically  visualized  as  2D  cross-sections  of  6D  stress  space.  The  “octahedral  diagram” 
in  Eigure  7,  for  example,  is  a  visualization  as  seen  looking  down  the  hydrostat  in  principal  stress 
space,  where  the  median  (50%)  quantile  surface  is  surrounded  by  other  quantile  surfaces  to  define  a 
“cloud”  that  replaces  the  discrete  limit  surface  in  deterministic  models.  The  “meridional  diagram” 
in  Figure  9  depicts  a  different  2D  slice  of  6D  stress  space  to  show  that  equivalent  shear  strength 
T  =  ^/72  increases  with  pressure  p  =  — /i/3,  where  /i  and  J2  are  the  first  and  second  invariants 
of  the  stress  tensor  T,  defined  Ii  —  trT  and  J2  —  |tr(devT)^.  Factors  of  ^/2  and  ^/3  in  the  axes 
labels  make  the  plot  isomorphic  to  stress  space  (ensuring  that  the  6D  unit  tensor  N  maps  to  a  2D  unit 
vector  that  also  forms  the  correct  angle  with  the  hydrostat).  In  upcoming  plots,  isomorphic  scaling 
factors  are  omitted  whenever  geometric  accuracy  is  not  needed.  Quantile  surfaces  in  Figure  9  are 
depicted  as  a  gray-scale  “cloud”.  The  red  fine  is  the  median  (50%  quantile),  taken  from  the  limit 
surface  in  a  deterministic  damage  model.  The  5%  and  95%  quantiles  are  lighter-gray  contours. 

The  loading  increment  direction  is  typically  fixed  in  model  calibration  experiments,  meaning  that 
the  stress  tensor  is  of  the  form  T  =  aM  -|-  C,  where  a  is  a  time  varying  load  intensity  parameter, 
M  is  a  unit  tensor  in  the  loading  increment  direction,  and  C  is  a  constant  tensor  that  is  usually  an 
elastic  prestress  reference  state.  If,  for  example,  the  experiment  involves  a  hydrostatic  compression 
to  a  pressure  p  followed  by  axisymmetric  compression  with  lateral  stress  held  constant  at  p,  then 


The  cloud  thickness 
(quantified  by  the  standard 
deviation)  may  vary  with 
ioading  obliquity. 


■dip 


Figure  9.  Loading  trajectory  through  a  fuzzy  failure  “cloud”. 


Copyright  ©  2014  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Engng  (2014) 
DOI:  10.1002/nme 


O.  E.  STRACK,  R.  B.  LEAVY  AND  R.  M.  BRANNON 


C  =  — DIAG[/>,  p,  p]  andM  =  DIAG[— 1, 0, 0].  If  C  =  0,  M  is  the  same  as  N  defined  in  Section  2.2 
and  sketched  in  Figure  9.  For  these  experiments,  strength  distributions  are  actually  distrihutions 
for  the  scalar  loading  parameter  a,  and  the  quantiles  for  a  are  then  simply  mapped  to  full  6D 
stress  space  to  generate  points  for  fitting  the  corresponding  quantile  surfaces.  If  a  is  the  median, 
for  example,  then  T  =  ctM  +  C  must  he  a  point  on  the  median  quantile  surface.  Points  on  other 
quantile  surfaces  in  6D  stress  space  are  similarly  identified.  The  process  of  finding  stress  tensors  on 
quantile  surfaces  must  be  repeated  for  at  least  as  many  distinct  loading  paths  as  there  are  distinct 
material  parameters  in  the  limit  function  F(T,  V),  thus  providing  a  means  of  setting  all  material 
parameters  corresponding  to  any  quantile  realization  of  strength.  The  strength  a  in  Equation  (7)  does 
not  necessarily  represent  a  single  component  of  stress  or  even  any  particular  stress  invariant.  It  is 
merely  a  scalar  parameter  in  the  range  from  0  to  oo,  whose  distribution  p(a,  V)  indirectly  quantifies 
variation  of  the  tensor-based  CDF  P^“’(T,  V)  along  a  fixed  loading  trajectory  M.  Specifically,  the 
generally  meaningless  scalar  a  is  mapped  to  the  physically  significant  stress  tensor,  T  =  aM  +  C. 
The  scalar-based  PDF  p{a,  V)  then  provides  information  about  the  tensor-based  CDF  V) 

through  the  directed  derivative  relationship,  p{a,  V)  —  (9P^“’/9T)  :  M. 

It  should  be  clear  from  examination  of  Figure  9  that  the  absolute  spread  in  failure  strength  data 
(quantified  by  variance  of  a)  will  be  affected  by  the  angle  of  attack  in  stress  space.  More  broadly, 
each  different  loading  direction  M  (i.e.,  each  test,  such  as  uniaxial  tension,  compression,  shear,  etc.) 
is  associated  with  a  different  distribution  of  strength.  They  are  connected  into  a  unified  theory  by  fit¬ 
ting  stress  tensors  (not  a)  on  each  quantile  surface  to  the  limit  function  F(T,  V)  of  a  “conventional” 
deterministic  damage  model. 

For  our  own  baseline  “conventional”  deterministic  damage  constitutive  theory,  we  created  a  new 
model,  called  Kayenta  [6],  which  revised  an  existing  geomaterial  model  [47]  to  incorporate  loss 
of  strength  in  a  manner  similar  to  the  Johnson-Holmquist  class  of  ceramic  models  [45,  48-50]. 
Kayenta  is  a  generalized  isotropic  damage  model  (as  opposed  to  an  anisotropic  model,  e.g.,  [28], 
or  a  combined  plasticity/linear-fracture-mechanics  model,  e.g.,  [51]).  Similar  to  the  strain-based 
Johnson-Cook  fracture  [46]  and  Johnson-Holmquist  models,  Kayenta  allows  equivalent  shear 
strength  t  =  in  axisymmetric  compression  to  depend  on  pressure  p  =  — /i/3  according  to 

T  =  ai  —  a3e~‘^^^  (8) 


where  {ai,  02,  <23}  are  empirical  fitting  parameters  (an  additional  term,  +04 p,  is  also  available  for 
linear  Drucker-Prager  (and  similar)  functions  but  is  omitted  here  for  simplicity). 

The  Johnson-Holmquist  and  Johnson-Cook  models  fall  into  the  category  of  nonlinear  Drucker- 
Prager  generalized  plasticity  (i.e.,  their  octahedral  profile  is  always  a  circle).  Kayenta  incorporates 
the  third  invariant  (the  Lode  angle  0)  by  multiplying  r  by  a  function  of  9  and  p  that  equals  unity  in 
axisymmetric  compression  but  is  larger  than  unity  at  other  Lode  angles  (giving  non-circular  octa- 
hederal  profiles  such  as  those  in  Figure  7).  In  Kayenta,  the  Lode  angle  multiplier  function  also 
includes  pressure  dependence  that  allows  the  octahedral  profile  to  vary  with  pressure,  smoothly  (and 
automatically)  transitioning  from  a  maximum  principal  stress  (or  strain)  criterion  at  low  pressures 
to  a  von  Mises  criterion  at  high  pressures.  A  cap  for  porous  media  is  also  available  as  an  option  but 
is  not  used  in  the  sample  calculations  herein. 

The  limit  surface  defined  by  Equation  (8)  could  be  perturbed  by  changing  one  or  more  of  the 
{a\,a2,  <23}  parameters.  However,  these  parameters  ought  not  to  be  varied  independently.  The  user 
must  ensure  that  PDFs  used  to  perturb  material  properties  will  not  violate  any  fundamental  physical 
constraints.  Moreover,  the  user  must  recognize  that  achieving  an  observed  Weibull  (or  other)  distri¬ 
bution  in  strength  might  require  complicated  non- Weibull  joint  distribution  functions  for  the  model 
parameters.  For  example,  the  constraint  ai  >  <23  must  be  satisfied  to  ensure  positive  shear  strength 
at  zero  pressure.  Other  considerations,  such  as  convexity  of  the  limit  surface,  place  additional  cou¬ 
pled  constraints  on  the  {ai,a2.‘23}  parameters.  None  of  these  strength-defining  parameters  are 
candidates  for  independent  Weibull  (or  similar)  perturbation.  To  address  this  inconvenience,  we 
re-parameterized  Equation  (8)  into  the  form 


r  =  T 


1  _  g-s(p+h)IY 


(9) 
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(a)  Possible  Meridional  view  (b)  Preferred  Meridional  view 


Figure  10.  Meridional  profiles  (including  probability  clouds  depicting  variability)  corresponding  to  (a) 
shear  dominant  variability  obtained  by  perturbing  Y ,  which  failed  to  match  observed  damage  patterns  in 
ceramic  fragmentation  simulations,  and  (b)  pressure-dominant  variability  obtained  by  perturbing  h,  which 
is  consistent  with  both  spall  and  Hugoniot  elastic  limit  data  and  gives  more  realistic  damage  patterns 

in  simulations. 

where  Y  is  the  shear  strength  limit,  h  is  the  spherical  tensile  strength,  and  s  is  the  initial  slope.  The 
meanings  of  these  alternative  parameters  are  labeled  in  Figure  10,  which  also  shows  probability 
density  clouds  resulting  from  statistically  perturbing  Y  or  h  individually. 

The  Weibull  modulus  for  SiC  was  reported  to  be  17  by  the  manufacturer  Coorstek- Vista  [52], 
but  spall  data  [53]  for  SiC-N  seem  to  have  a  Weibull  modulus  of  only  4  (i.e.,  much  more  variable). 
Brazilian  tensile  strengths  (presented  later  in  Figure  17)  appear  to  have  a  Weibull  modulus  of  about 
8.  Flugoniot  data  for  B4C  [54]  has  a  Weibull  modulus  of  60,  with  similar  numbers  for  SiC.  Moreover, 
it  is  not  even  clear  that  some  of  these  data  sets  are  very  well  fitted  to  a  Weibull  distribution.  To  match 
the  general  character  of  observed  strength  data  (i.e.,  lessening  of  variability  with  pressure  and  non- 
Weibull  strength  distributions  for  most  loading  trajectories),  it  is  proposed  that  the  limit  surface  be 
perturbed  as  sketched  in  Figure  10(b). 

Unlike  the  original  {ai,  0^,03}  fitting  parameters,  there  is  no  coupling  of  the  constraints  on 
Kayenta’s  alternative  limit-surface  parameters  {T,  h,  j},  which  only  must  be  positive.  Flence,  any  of 
them  may  be  Weibull  perturbed  independently  according  to  Equation  (7)  (where  a  in  that  equation 
is  replaced  with  Y,  h,  or  s).  Consistent  with  a  similar  statistical  uncertainty  assumption  used  to  ana¬ 
lyze  the  nearly  spherical  tension  states  that  exist  during  spall  [55],  and  also  consistent  with  the 
tendency  (evident  in  Figure  13)  of  strength  measurements  to  be  more  variable  at  low  pressure,  we 
take  the  h  parameter  to  be  Weibull  distributed,  whereas  the  other  limit  surface  parameters  (Y  and 
s)  are  deterministic.  This  choice  leads  to  non- Weibull  statistics  and  non- Weibull  scale  effects  for 
all  non-spherical  loading  paths.  This  approach,  which  truly  accounts  for  the  6D  tensorial  nature  of 
stress,  differs  significantly  from  related  efforts  [56]  that  apply  a  ID  Weibull  distribution  to  a  single 
stress  invariant  such  as  equivalent  shear  stress,  maximum  principal  stress,  and  so  on.  We  conjecture 
that  variability  decreases  with  pressure  because  the  additional  confinement  provides  friction  that 
can  lock  crack  faces  and  hence  suppress  brittle  failure  modes  in  shear.  Of  course,  high-pressure  data 
often  comes  from  high-rate  experiments,  but  we  have  tentatively  ruled  out  the  possibility  that  vari¬ 
ability  decreases  with  increasing  loading  rate;  this  choice  is  motivated  by  the  fact  that  quasistatic 
indentation  experiments  also  seem  to  exhibit  pressure-based  variability  because  they  tend  not  to 
bifurcate  into  symmetry-breaking  radial  cracks  in  regions  of  high  confinement  such  as  immediately 
beneath  an  indenter  [57].  Further  insight  into  the  influence  of  pressure  and  loading  rate  is  provided 
in  a  recent  investigation  of  rate  effects  in  brittle  media  [58],  which  applies  micromorphological 
dimensional  analysis  to  obtain  a  universal  curve  which  closely  fits  rate-dependent  and  statistically 
variable  strength  for  several  different  cracked  and  porous  brittle  materials  (including  various  rocks, 
ceramics,  and  concrete). 

In  Figure  1 1(a),  the  value  of  the  spherical  tension  limit  h  is  exaggerated  to  illustrate  the  character 
of  perturbation  effects.  The  actual  value  of  h  for  brittle  ceramics  is  typically  so  small  in  comparison 
with  the  high  pressure  strength  Y  that  the  actual  median  limit  curve  might  appear  to  pass  through 
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Figure  11.  The  strength  size  effect  in  “directed- WeibuH”  theory  varies  with  stress  space  trajectory.  Here, 
strength  can  increase  without  bound  for  progressively  smaller  samples  in  spherical  tension  (horizontal 
green  loading  line  in  (a)  and  straight  green  line  in  (b)),  but  it  naturally  saturates  to  a  classical  plasticity 
limit  in  uniaxial  tension  (black  angled  loading  line  in  (a)  and  black  curved  line  in  (b)),  with  clear 
deviation  from  Equation  (3)  for  non-spherical  loading.  This  general  framework  naturally  accommodates 
any  other  experimentally  observed  distributions  of  quantile  surfaces,  including  non-monotonic  scaling 

of  strength. 


the  origin  for  lab  scale  (cm^)  samples,  as  illustrated  in  Figure  13(a)  for  a  large  (>  1  cm^)  sample 
of  SiC-N.  Using  Figure  7  to  generate  strength  envelopes  means  that  realizations  of  the  {Y,h,s} 
parameters  are  scale  dependent.  Smaller  samples  have,  on  average,  higher  strength.  To  illustrate. 
Figure  13(h)  shows  the  median  limit  surface  for  a  grain-scale  (/xm^)  sample. 

When  only  the  h  parameter  is  Weihull-perturhed,  the  median  spherical  tensile  strength  approaches 
infinity  as  the  element  size  approaches  zero.  However  (as  seen  in  Figure  11),  the  median  uniaxial 
tension  and  shear  strengths  have  finite  limits,  and  the  strength  realizations  in  these  non-spherical 
directions  deviate  realistically  from  a  Weihull  distribution.  This  “directionally-Weibull”  approach 
gives  a  scaling  effect  similar  to  that  reported  by  Bazant  [23]  on  the  basis  of  experimental  data 
for  concrete. 

Bazant  has  amassed  direct  experimental  data  for  concrete  supporting  deviation  from  the  Weihull 
size  dependence  similar  to  what  is  depicted  in  Figure  11(b).  No  similar  direct  evidence  of  the 
deviation  from  the  Weihull  size  effect  has  (to  our  knowledge)  been  collected  for  sufficiently  small 
specimens  of  structural  ceramics  (whose  “aggregate”  size  is  orders  of  magnitude  smaller  than  that 
of  concrete).  We  arrived  at  this  effect  indirectly  by  recognizing  that  macroscale  tensile  strengths  are 
far  more  variable  than  compressive  strengths  (Figure  13),  suggesting  that  limit  surface  realizations 
require  perturbations  of  the  type  in  Figure  10(b). 

The  nature  of  the  limit  surface  “variability  cloud”  and  its  median  master  surface  strongly  influ¬ 
ences  the  appearance  of  cone  or  radial  cracks  in  dynamic  indentation.  In  indentation  experiments 
(discussed  in  Section  8),  cone  cracks  tend  to  emerge  around  and  under  the  area  of  contact  where 
the  pressure  is  extremely  high.  They  are  also  symmetric  [59],  indicating  a  low  degree  of  spatial 
variability  in  the  high-pressure  realm.  Radial  cracks,  on  the  other  hand,  tend  to  appear  in  rela¬ 
tively  low-pressure  regions,  and  they  significantly  disrupt  symmetry,  indicating  a  high  degree  of 
spatial  variability  in  strength  akin  to  Figure  10(b).  As  seen  in  Figure  12,  similar  observations  have 
been  made  for  brittle  coatings  [57,  60],  glass  [61],  and  boron  carbide  plates  [62].  The  order  of  the 
appearance  of  radial  cracks  and  cone  cracks  can  be  altered  by  changes  in  materials  and  indenter 
geometries. 

Mapping  out  aleatory  quantile  surfaces  is  clearly  labor  intensive.  Given  that  it  is  often  impossible 
or  impractical  to  obtain  statistically  significant  samples  of  laboratory-measured  strengths,  judicious 
use  of  virtual  data  (e.g.,  [64,  65])  might  help  to  complete  the  aleatory  model  parameterization  and/or 
to  identify  appropriate  functional  forms  of  the  limit  function  F(T,  V)  about  which  perturbations 
must  be  applied.  Virtual  data  (if  validated  against  available  data  and  verified  against  classical  ideal¬ 
ized  exact  solutions)  provides  insight  into  the  microphysical  sources  of  aleatory  uncertainty,  which 
likewise  aids  in  the  formulation  of  path-dependent  macroscale  models  for  evolution  of  internal  state 
variables  that  control  the  loss  of  stiffness  and  collapse  of  the  limit  surface  from  failure. 
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Figure  12.  (a)  Simulated  dynamic  indentation:  transition  from  conical  damage  to  radial  cracking  beneath 
the  cone,  similar  in  character  to  observations  in  glass  indentation  [57,  60].  (b)  Quasi-static  indentation: 
symmetry-breaking  radial  cracks  beneath  symmetry-preserving  ring  cracks  [61,  63]. 

5.  MODEL  CALIBRATION  OF  VARIABILITY  AND  SCALE  EFFECTS 


Methods  to  account  for  strength  data  from  differently  sized  experimental  samples  and  methods  to  fit 
repeated  strength  measurements  to  a  Weibull  (or  similar  size-dependent)  distribution  are  discussed 
in  this  section.  These  techniques  provide  median  strength  values  and  distribution  properties  (such 
as  the  mean  and  standard  deviation)  for  a  fixed  experimental  configuration. 

5.7.  Normalizing  experimental  data  to  a  standard  reference  volume 

Ideally,  the  same  size  sample  would  be  used  in  a  variety  of  strength  tests  for  various  loading  modes 
such  as  triaxial  compression,  simple  shear,  buttonhead  tension,  and  so  on.  If  median  strengths  are 
measured  in  these  different  loading  modes  using  the  same  sample  size,  then  the  results  directly  map 
to  a  strength  limit  surface  in  stress  space.  However,  if  different  loading  modes  are  measured  using 
samples  of  different  sizes,  then  all  of  the  data  must  be  scaled  to  a  common  reference  volume  before 
mapping  the  median  strength  surface. 

Weibull  scaling  in  Equation  (2),  including  the  reformulated  versions  in  Equation  (7),  and 
Equation  (11),  requires  the  volume  to  be  expressed  relative  to  a  reference  volume  V.  Similarly,  the 
median  stress  a  corresponds  to  measurements  using  a  sample  of  size  V.  Recognizing  that  different 
experiments  might  be  performed  using  samples  of  different  volumes,  this  section  briefly  explains 
how  all  data  from  homogeneous  loading  experiments  can  be  scaled  to  a  common  reference  volume. 
The  far  more  subtle  problem  of  scaling  data  from  heterogeneous  loading  experiments  (such  as 
Brazilian  or  bending  tests)  remains  unsolved,  but  the  relevant  issues  are  introduced  in  Section  6. 

Consider  a  set  of  experiments  that  measure  a  median  strength  a*  and  Weibull  modulus  m*  for 
a  sample  volume  V* .  According  to  Equation  (2),  the  same  set  of  experiments  conducted  using  a 
sample  with  the  reference  volume  V  would  produce  the  same  Weibull  modulus  (m  —  m*)  and  a 
volume-scaled  median  strength  given  by 


By  using  this  formula,  experiments  that  use  different  specimen  sizes  may  all  be  standardized  to  a 
single  reference  volume  V .  The  choice  for  the  reference  volume  is  arbitrary,  but  a  reference  vol¬ 
ume  is  essential  to  generate  failure  threshold  realizations  by  comparing  the  finite-element  size  to 
the  reference  volume  size  used  in  laboratory  calibration  experiments.  In  geomechanics,  it  is  com¬ 
mon  to  employ  two  different  sets  of  material  properties,  one  for  the  so-called  laboratory  scale  and 
another  for  “held  scale”  (cf.  [66]).  The  outlined  implementation  of  Weibull  theory  in  Equation  (10) 
automates  scale  effects  so  that  only  one  material  data  set  spans  a  large  range  of  size  scales. 

The  directionally- Weibull  theory  developed  in  Section  4  produces  a  different  apparent  (and 
approximate)  Weibull  modulus  m  for  different  loading  trajectories,  which  should  be  accounted  for 
when  scaling  sample  sizes.  The  Weibull  modulus  for  a  selected  model  and  loading  direction  can  be 
determined  by  numerically  (or  analytically)  generating  a  series  of  realizations  of  the  limit  surface. 


(10) 
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collecting  the  corresponding  CDF  (failure  probability)  values  along  the  fixed  loading  trajectory,  and 
then  applying  the  analysis  outlined  below. 

5.2.  Measuring  strength  parameters 

The  following  procedure  for  measuring  the  strength  parameters  is  performed  under  two  idealized 
conditions.  One  is  that  the  loading  is  statistically  homogeneous  with  no  significant  macroscale  stress 
gradients  prior  to  failure  (e.g.,  quasi-static  triaxial  compression  and  extension  and  buttonhead  ten¬ 
sion  tests).  The  second  condition  presumes  that  the  same  reference  specimen  size  V  is  used  in  all 
experiments.  Even  though  these  conditions  can  be  difficult  to  meet,  understanding  the  method  used 
to  calibrate  data  under  these  idealizations  is  a  prerequisite  to  subsequent  discussions  in  Section  6 
that  release  these  assumptions  by  accounting  for  stress  gradients  within  the  specimen. 

Although  expanded  testing  and  analysis  are  likely  to  be  ongoing  for  many  years,  preliminary 
scale-normalized  strength  data  for  SiC-N  are  shown  in  Figure  13  for  a  variety  of  pressures,  with 
median  scaled  data  fitted  to  Kayenta’s  limit  surface  function.  Without  scaling  all  data  to  a  common 
specimen  size,  no  fitting  exercise  would  be  meaningful.  The  raw  data  measured  using  a  sample  of 
volume  V  were  scaled  to  a  common  reference  volume  V  via  Equation  (10).  Eor  experiments  having 
a  nonuniform  stress  field,  the  specimen  volume  V  is  understood  to  be  the  effective  volume  from  a 
data  delocalization  scheme  (Section  6). 

A  suite  of  experiments  using  different  loading  directions  must  be  run  to  identify  the  median 
failure  envelope,  as  it  is  loading  path-dependent.  This  conventional  part  of  the  calibration  process 
identifies  the  dependence  of  a  on  stress  triaxiality  (or,  more  generally,  on  N).  Once  the  loading  path 
for  a  particular  experiment  is  selected,  the  experiment  must  be  conducted  multiple  times  to  obtain 
a  statistically  significant  collection  of  strength  realizations  {ai ,  02, . . . ,  CTm}.  Deciding  if  the  data  fit 
a  Weibull  distribution  requires  assigning  values  to  the  reference  volume  V ,  the  median  strength  a 
at  the  reference  volume  V,  and  the  Weibull  modulus  m.  Specifically,  the  data  must  be  fit  to  the 
“linearized”  version  of  the  Weibull  equation 

y  —  mx  +  b  where  v  =  lnln( - lnln(2),  x  —  and  6  =  ln=,  (11) 

ypsatej  (j  y 

which  is  obtained  by  taking  the  natural  logarithm  of  Equation  (2)  twice  and  rearranging.  Any  data, 
whether  Weibull  distributed  or  not,  may  be  depicted  in  a  plot  of  y  versus  x.  Such  data  are  proved  to 
be  Weibull  distributed  if  the  data  fall  on  a  straight  line  with  respect  to  these  Weibull  axes,  though  we 
caution  that  even  exactly  Weibull-distributed  data  will  exhibit  stair-stepped  deviations  from  the  line 
as  a  result  of  finite-sampling  errors  (Eigures  14  and  16).  As  seen  in  Eigure  16,  finite  sampling  errors 
in  exactly  Weibull-distributed  data  tend  to  produce  a  “dip”  away  from  the  lower-left  portion  of  the 
theoretical  Weibull  line.  On  the  other  hand,  simple  microphysically  derived  strength  CDEs  (cf.  [31]) 
also  have  a  dip  that  is  not  the  result  of  finite  sampling  error.  Resolving  whether  a  dip  is  a  genuine 


(a)  Data  clouds  and  experimental  data  sealed  to  a  large  (b)  Data  elouds  and  notional  experimental  data  scaled  to 
sample  reference  volume  V.  a  small  sample  reference  volume  V. 


Figure  13.  Examples  of  strength  scaling  and  data  clouds  for  comparison  with  experiments  of  silicon  carbide, 
a)  Data  clouds  and  experimental  data  scaled  to  a  large  sample  reference  volume  V ;  (b)  data  clouds  and 
notional  experimental  data  scaled  to  a  small  sample  reference  volume  V . 
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Figure  14.  Verification  of  mesh  independence  for  the  onset  of  failure.  The  figures  show  simulated  (jc*,  yk) 
data  from  Equation  (13)  for  three  mesh  resolutions  to  demonstrate  that  failure  statistics  for  the  onset  of 
failure  are  mesh  insensitive  and  approach  the  Weibull  distribution  as  the  number  of  realizations  is  increased. 

(a)  600  points  and  (b)  6000  points. 


attribute  of  the  failure  statistics  or  merely  a  byproduct  of  finite  sampling  is  a  strong  motivation  to 
increase  the  number  of  strength  tests  in  laboratory  testing. 

As  discussed  in  Sections  2.2  and  4,  “strength”  is  regarded  as  the  value  of  the  loading  parameter 
a  along  a  loading  trajectory  at  failure.  Here,  “failure”  is  the  initiation  of  stiffness  and  strength 
reduction  upon  attainment  of  a  peak  value  of  a,  but  the  definition  could  be  altered  to  mean  phase 
transformation,  or  any  other  “event”  marking  a  distinct  irreversible  change  in  material  behavior. 
Repeating  the  experiment  n  times  will  give  a  set  of  strength  realizations  that  can  be  put  into  canon¬ 
ical  order,  ai  ^  ai  ^  . . .  ff*  •  •  •  ^  ffn  ■  Associated  cumulative  probabilities  for  sample  number  k  are 
given  by 

pfeii  ^  ^  ^ 

This  estimator  is  advocated  by  Sullivan  and  Lauzon  [67],  who  also  recommend  using  no  fewer  than 
20  measurements  and  preferably  more  than  50.  An  even  larger  number  of  measurements  is  required 
to  determine  if  the  Weibull  distribution  is,  in  fact,  the  actual  distribution.  For  example,  more  than 
1000  tests  would  be  required  to  distinguish  between  the  Weibull  distribution  and  the  log-normal 
distribution  [68].  Other  possibilities  such  as  the  Gumbel  distribution  [69]  undoubtedly  require  an 
equally  large  sample  to  validate.  Consequently,  no  part  of  this  paper  presupposes  that  a  Weibull 
distribution  applies  to  measured  strength  distributions.  A  Weibull  distribution  is  used  primarily  as 
a  familiar  example  to  illustrate  considerations  in  data  fitting.  As  explained  in  Section  4,  the  empir¬ 
ical  framework  naturally  predicts  strengths  fitting  no  standard  distribution,  Weibull  or  otherwise. 
Minimally,  distributions  must  be  used  that  will  generate  physically  admissible  realizations. 

For  fitting  data  to  a  Weibull  distribution,  the  volumes  V  and  V  are  both  set  equal  to  the  specimen 
volume  used  in  the  experiment,  making  6  =  0  in  Equation  (11).  From  the  collection  of  (a^, 
data  pairs,  a  is  identified  to  be  the  median  (i.e.,  the  value  corresponding  to  =  1  /2).  Thus,  only 
the  Weibull  modulus  m  in  Equation  (11)  remains  yet  to  be  inferred  from  the  measured  data.  To  obtain 
m,  the  (Uyt,  data  pairs  are  converted  to  {Xk,yk)  data  pairs  by  the  following  transformations: 

iTt  1 

Xyt  =  In  —  and  =  Inin - 5 — In  In 2.  (13) 

a  F| 

The  Weibull  modulus  m  is  then  determined  via  true  linear  (not  affine)  regression  fitting  of  these 
data  pairs  to  y  =  mx.  Namely, 

n 

yn^n 

m  =  ^ -  (14) 

k=l 
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Naturally,  the  regression  error  is  a  direct  quantitative  indicator  of  the  extent  to  which  a  Weihull 
distribution  is  appropriate.  Sullivan  and  Lauzon  [67]  offer  an  alternative  to  simple  linear  regression. 
However,  hy  analyzing  synthetically  generated  exactly  Weibull-distrihuted  data,  we  find  their 
method  to  be  superior  to  linear  regression  only  when  the  total  number  of  strength  measurements /ar 
exceeds  their  recommended  minimum  of  20  data  points.  For  data  sets  of  fewer  than  14  points,  we 
similarly  find  that  Equation  (14)  more  accurately  predicts  the  Weihull  modulus  than  a  scheme  that 
sets  the  two  (shape  and  scale)  Weihull  parameters  by  matching  the  sample  mean  and  standard  devi¬ 
ation  to  Equations  5  and  6.  Even  more  than  14  points  are  needed  for  this  alternative  fitting  scheme 
to  be  superior  whenever  data  beyond  2  or  3  standard  deviations  from  the  mean  are  misidentified 
as  “outliers”.  Eor  a  typical  engineering  sample,  which  rarely  exceeds  20  data  points  and  is  often 
far  below  what  is  needed  for  statistical  significance,  these  finite-sampling  errors  and  inadvertent 
removal  of  outliers  seem  to  be  best  managed  using  the  standard  CDE  estimator  in  Equation  (12)  and 
the  simple  linear  regression  in  Equation  (14).  Overall,  the  median  is  the  most  convenient  measure 
of  central  tendency  because  our  approach  aims  to  identify  generally  non- Weihull  quantile  surfaces 
(for  which  the  median  is  the  50%  quantile). 

Eigure  14  verifies  mesh  independence  (and  recovery  of  the  user-prescribed  Weihull  strength  dis¬ 
tribution)  for  the  onset  of  failure  under  idealized  spherical  tension.  The  effect  of  finite  sampling  on  a 
Weihull  plot  (i.e.,  a  plot  of  j  vs.  x)  is  apparent  in  the  numerical  data;  finite  sampling  tends  to  cause 
dips  away  from  the  Weihull  line  near  the  tail  of  the  distribution  (see  also  Eigure  16). 

In  Eigure  15(a),  a  ceramic  flyer  (the  darker  color  on  right  side  of  both  images)  impacts  a  ceramic 
target  twice  the  thickness  of  the  flyer.  The  damaged  (spalled)  material  has  been  hidden  from  the 
visualization  of  the  simulation  on  the  left,  and  the  velocities  are  shown  on  the  right  (blue  to  red 
color  scale  indicates  low  to  high  velocities,  respectively).  Spatial  variability  induces  non-uniform 
velocity  on  the  surface  of  the  target  that  would  not  occur  in  a  deterministic  simulation.  Spall  is 
a  dynamic  tensile  failure  produced  by  the  interaction  of  two  rarefaction  release  waves,  making  it 
a  truly  internal  failure  mode  that  is  unaffected  by  surface  smoothness  at  the  specimen  boundary. 
Because  of  the  strongly  non-uniform  stress  fields  occurring  in  spall  tests,  the  data  delocalization 
methods  discussed  in  Section  6  must  be  applied  to  iterate  on  the  specimen  volume  associated  with 
measured  spall  strengths. 

In  Eigure  16,  Weihull  plots  and  plots  of  the  complementary  CDF  are  shown  for  two  different 
experiments:  the  Hugoniot  elastic  limit  strengths  for  B4C  ceramic  [54]  and  previously  unpublished 
tensile  strengths  for  SiC-N  for  the  Brazilian  diametrical  compression  test.  Both  appear  to  be  Weibull- 
distrihuted,  but  the  B4C  shock  compression  data  clearly  exhibit  much  lower  variability  than  the 
Brazilian  data.  These  data  deviate  from  the  best  Weihull  fit  but  only  by  amounts  well  within  the 
range  of  deviations  of  Weihull  data  pulled  from  the  exact  analytical  Weihull  distribution  for  the  same 
modulus  m  (also  shown  in  the  figure).  As  seen,  finite-sampling  deviations  in  the  exactly  Weibull- 
distrihuted  data  are  similar  in  magnitude  and  character  to  deviations  of  the  measured  data  from  the 
Weihull  fit.  Thus,  measured  data  may  be  closer  to  a  Weihull  distribution  than  might  seem  evident 
from  a  Weihull  plot. 


(b)  Brazillian  test 


(a)  Spall  test 


Figure  15.  Simulations  with  variability  demonstrating  bifurcation  from  nominal  symmetries,  (a)  Spall  test 

and  (b)  Brazilian  test. 
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psafe 


Figure  16.  Measured  strength  data  (thick  black  line)  and  best  Weibull  fit  (thick  straight  gray  line),  along 
with  the  same  number  of  exactly  Weibull  distributed  simulated  data  points  (thin  colored  lines). 


ln(o-^“')  ln(a^“') 

(a)  Brazilian  experiments  (b)  Brazilian  simulations 

Figure  17.  (a)  Laboratory  strength  data  show  effects  of  three  specimen  sizes  (five  strengths  measured  for 
the  largest  specimen  are  near  the  left-most  linear  regression  fit,  twenty  strengths  for  the  medium-sized 
specimen  are  in  the  center,  and  six  for  the  smallest  specimen  are  on  the  right),  (b)  Simulations  for  the 
same  sample  sizes,  with  initial  spherical  strength  perturbed  using  a  Weibull  modulus  of  6,  update  those 

in  [70]. 


Finite-sampling  effects  (the  “dipping”  deviation  from  the  Weibull  line)  are  especially  apparent 
in  the  Brazilian  strength  data  plotted  in  Figure  17,  because  those  experiments  acquired  fewer  data 
points  than  the  experiments  in  the  left  side  Flugoniot  elastic  limit  data  of  Figure  16.  Examples  of 
Brazilian  experiments  and  comparable  simulations  are  in  Figure  15(b). 

The  strength  size  scaling  of  the  Brazilian  material  samples  is  noted  in  both  the  experiments  and 
simulations  in  Figure  17.  These  plots  omit  the  typical  normalization  of  the  abscissa  by  the  median 
stress  at  failure  a  so  that  the  relative  strengths  and  slopes  can  more  readily  be  seen  for  the  different 
sample  sizes.  An  interesting  aspect  of  Figure  17  is  that  the  experimental  data  appear  to  show  not  only 
an  increase  in  strength  with  decreasing  specimen  size  but  also  a  decreasing  Weibull  modulus  (i.e., 
decreasing  slope  and  hence  greater  relative  variability  in  strength)  with  decreasing  specimen  size. 
However,  because  of  the  small  number  of  experimental  samples  and  the  impact  of  finite-sampling 
errors,  this  trend  is  inconclusive.  The  numerical  simulations  in  Figure  17  do  not  show  the  same 
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trend,  and  further  work  is  warranted  to  understand  the  scaling  of  the  apparent  Weihull  modulus  as  a 
function  of  sample  size  for  the  Brazilian  test  (both  experimentally  and  numerically). 

Although  it  would  he  possible  to  tune  the  parameters  of  the  model  to  better  match  the  data,  doing 
so  would  be  grossly  misleading,  because  the  simulations  failed  to  exhibit  the  same  mitigation  of 
mesh  sensitivity  that  was  seen  in  the  dynamic  indentation  simulations  of  Figure  1  (i.e.,  changing 
the  mesh  would  require  inadmissible  re-tuning  of  material  properties).  In  the  Brazilian  simulations, 
the  methods  of  Section  6  have  been  confirmed  to  provide  a  high  rate  of  convergence  for  the  onset 
of  failure  in  the  Brazilian  test.  However,  further  research  breakthroughs  are  required  to  likewise 
identify  regularizations  for  the  subsequent  progression  of  failure  (both  at  the  element  and  structural 
level)  while  also  predicting  realistic  strength  distributions  and  scale  effects. 

6.  NONLOCAL  CORRECTIONS  OF  DISCRETIZATION  ERRORS 

It  has  long  been  asserted  that  finite-element  analysis  of  distributed  damage  requires  some  type  of 
nonlocal  theory  [71,  72].  Nonlocality  is  attributed  to  many  physical  sources  such  as  crack  inter¬ 
actions  [73],  which  might  be  important  on  fine  discretizations  where  aleatory  uncertainty  in  crack 
locations  implies  the  need  for  correlated  strength  distributions.  However,  recent  work  [32]  suggests 
purely  numerical  motivations  for  nonlocality. 

When  attempting  to  calibrate  scale-dependent  strength  from  laboratory  experiments,  so-called 
data  delocalization  techniques  are  recommended  to  assign  the  corresponding  specimen  volume.  If 
the  stress  field  is  uniform,  then  the  specimen  volume  is  the  actual  volume.  However,  for  bend¬ 
ing,  Brazilian  indirect  tension,  spall,  and  other  experiments  having  a  nonuniform  stress  field, 
data-delocalization  theory  [38,  74-76]  sets  the  specimen  volume  using  a  stress-weighted  aver¬ 
age.  For  a  uniformly  loaded  and  statistically  homogeneous  specimen,  we  adopt  the  axiom  that 
the  probability  of  a  local  failure  must  scale  with  volume,  not  area  (i.e.,  with  V,  not  E^/^).  How¬ 
ever,  our  own  experimental  data  for  the  Brazilian  indirect  tension  tests  in  Figure  17  scale  with 
area,  not  volume.  This  apparent  area-based  scaling  is  a  natural  result  when  the  data  are  analyzed 
using  data-delocalization  theory.  Specifically,  whenever  stress  is  nonuniform  (so  that  only  a  por¬ 
tion  of  the  actual  volume  V  is  loaded  to  an  appreciable  stress  level),  the  “effective  volume”  is 
defined  by 


Eeff=  f  ia/ap,,0'"dV  (15) 

Jv 

where  a  is  a  scalar  measure  of  stress  (such  as  the  maximum  principal  stress)  and  Opeak  is  the 
maximum  of  a  in  the  domain  V.  In  the  Brazilian  test,  for  example,  higher  stresses  are  located  over  a 
narrow  region  at  the  center  of  the  domain,  making  this  region  have  higher  weight  in  Equation  (15), 
thus  predicting  that  Egff  scales  with  area  as  observed  in  the  experiments. 

Of  course,  evaluating  Eeff  requires  a  description  of  the  non-uniform  stress  field  prior  to  failure. 
Iterations  are  also  required  to  find  Eeff  because  it  depends  on  one  of  the  properties  being  measured 
(the  Weihull  modulus,  m).  The  explicit  presence  of  the  Weihull  modulus  corresponds  to  an  implicit 
assumption  that  the  strength  is  Weihull  distributed.  Furthermore,  because  Equation  (15)  depends  on 
a  single  scalar  stress  invariant,  it  fails  to  account  for  the  tensorial  nature  of  stress,  and  it  implicitly 
introduces  a  failure  criterion  into  a  situation  where  the  failure  criterion  is  part  of  the  information 
being  sought  in  the  experiment  itself.  To  alleviate  some  of  these  shortcomings,  we  expect  that  a 
practical  generalization  of  Equation  (15)  might  be 

^eff=  f  (16) 

Jv  ^(A)peak 

where  the  form  and  function  of  X(T)  should  be  regarded  as  the  primary  information  sought  in  lab¬ 
oratory  testing.  As  explained  in  an  idealized  context,  namely  Equation  (11)  of  [31],  the  generalized 
function  “X”  is  expected  to  be  a  function  of  the  loading  mode  and  a  functional  of  underlying 
micromorphology  (e.g.,  evolving  flaw  size  and  orientation  statistics). 
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Data-delocalization  is  required  to  interpret  data  from  experiments  having  a  nonuniform  stress 
field  to  assign  a  volume-dependent  strength  under  uniform  stress.  Kamojjala  and  Brannon  [32]  have 
developed  a  “data-relocalization”  technique  to  reduce  errors  due  to  under-resolution  of  stress  fields. 
They  explain  that,  to  determine  the  failure  initiation  probability,  the  effective  uniform  stress  acting 
on  a  finite  element  should  be  set  as  an  m-norm  of  the  actual  pre-failure  stress  field: 


safe  _ 


where 


(Teff  = 


///a-dF 


1/m 


(17) 


Here,  m  is  the  Weibull  modulus,  implying  that  this  theory  would  require  revision  for  non- 
Weibull  distributions,  perhaps  replacing  a"®  with  a  micro-physically  based  X(T)  function,  as  in 
the  delocalization  proposal  of  Equation  (16).  Also,  a  is  the  pre-failure  analytical  principal  stress 
invariant  field  over  the  finite  element,  which  is  known  for  the  Brazilian  test  and  the  other  simple 
cases  studied  by  Kamojjala  and  Brannon.  To  adapt  this  relocalization  formula  to  cases  for  which  the 
analytical  stress  field  is  unknown,  a  nonlocality  algorithm  is  required  to  estimate  the  stress  concen¬ 
trations  within  an  element  on  the  basis  of  stresses  in  surrounding  elements.  This  recommendation  is 
a  purely  numerical  motivation  for  nonlocality,  and  it  might  be  accommodated  in  a  “semi-local”  way 
(i.e.,  without  needing  to  query  fields  in  neighboring  elements)  by  applying  an  m-norm  to  the  non¬ 
constant  stress  field  available  within  higher-order  elements  (Kamojjala  and  Brannon’s  exploration 
of  higher-order  elements  essentially  missed  this  opportunity  by  treating  the  stress  field  as  constant 
in  the  neighborhood  of  the  Gauss  points). 

The  previous  sections  have  focused  on  statistics  for  the  onset  of  failure,  which  controls  the 
statistically  variable  strength  assigned  at  initialization  and  data-relocalization  methods  to  adjust 
the  strength  of  an  element  in  response  to  stress  gradients.  In  a  rigorous  analysis  of  quasistatic 
mode-I  opening  of  a  sharp  crack,  Bazant  [23]  addresses  additional  considerations  associated  with 
failure  progression  that  may  or  may  not  lead  to  cascading  structural  failure  after  failure  initiation. 
Specifically,  Bazant  pointed  out  that  stress  singularities  near  a  crack  tip  increase  more  rapidly  with 
mesh  refinement  than  element  strengthening  associated  with  the  size  effect,  thus  encouraging  crack 
growth.  The  related  topic  of  damage  progression  is  the  subject  of  the  next  section. 


7.  STRENGTH  AND  STIEFNESS  DEGRADATION 

The  preceding  sections  focused  almost  exclusively  on  the  need  to  assign  different  values  of  material 
strength  to  each  element  in  a  computational  domain.  The  strength  assigned  to  different  elements 
is  statistically  variable  and  size-scale  dependent,  consistent  with  laboratory  observations.  This 
variability  in  strength  provides  perturbations  in  the  stress  field  necessary  to  stimulate  inherent 
instabilities  such  as  bifurcation  of  axisymmetric  loading  into  radial  cracking.  However,  once  cas¬ 
cading  failure  has  begun,  additional  scale  effects  are  well  known  to  be  necessary  to  mitigate  mesh 
dependence  during  the  progression  of  failure  [77]. 

Seminal  cohesive  zone  formulations  [78,  79]  and  more  recent  models  (cf.  [30,  80,  81])  introduce 
a  scale  effect  by  requiring  a  certain  amount  of  energy  be  required  to  advance  a  crack.  Until  recently, 
cohesive  zone  formulations  have  been  developed  on  the  basis  of  an  implicit  assumption  that  a  finite 
domain  contains  exactly  one  localized  zone.  However,  high-rate  loading  can  produce  a  network  of 
actively  growing  cracks  within  a  domain,  therefore  calling  for  an  increase  in  required  fracture  energy 
[58,  77].  As  explained  in  this  section,  this  overall  effect  is  accommodated  in  our  model  by  requiring 
a  finite  time  to  elapse  before  an  element  can  lose  strength. 

To  simulate  crack  growth,  conventional  isotropic  damage  models  usually  employ  a  scalar  damage 
parameter,  D,  that  varies  from  0  for  undamaged  material  to  1  for  fully  failed  material  [44].  These 
models  lead  to  a  drop  in  both  stiffness  and  strength  with  increasing  strain  (i.e.,  softening)  and 
accordingly  induce  localization  of  material  response.  Without  such  localization,  the  governing 
equations  can  change  type,  potentially  requiring  solution  of  a  diffusion  equation  [82].  On  the  other 
hand,  traditional  viscous  rate  dependence  in  the  governing  equations  can  have  a  regularizing  effect 
[83].  Similarly,  the  extreme  localization  of  brittle  and  quasi-brittle  materials  (i.e.,  formation  of 
macrocracks)  allows  the  governing  equations  to  remain  hyperbolic  as  long  as  boundary  conditions 


Copyright  ©  2014  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Engng  (2014) 
DOI:  10.1002/nme 


O.  E.  STRACK,  R.  B.  LEAVY  AND  R.  M.  BRANNON 


at  the  new  surfaces  are  handled  properly.  In  a  smeared  damage  model,  these  boundary  conditions 
are  approximated  hy  a  boundary  layer  of  failed  material  that  offers  resistance  to  compression  but 
progressively  reduced  resistance  to  tension  and  ultimately  only  frictional  resistance  to  shear.  Fur¬ 
ther  discourse  about  the  character  of  conventional  damage  models  is  beyond  the  scope  of  this  paper, 
where  our  focus  here  is  on  seeding  these  models  realistically  with  uncertainty  and  scale  effects.  The 
conventional  damage  model  used  for  our  case  study  (Kayenta)  is  well  documented  [6],  so  we  merely 
call  attention  to  some  aspects  that  differ  from  similar  models  available  in  commercial  codes. 

The  Kayenta  damage  model  employs  a  distinctive  approach  to  degrading  strength.  Taking  a  cue 
from  experimental  observations  [84]  that  cracks  in  high-stress  environments  grow  at  approximately 
constant  speed  (the  effect  of  which  is  seen  in  Figure  18),  the  model  smoothly  decreases  strength 
and  elastic  stiffness  to  failed  values  over  a  scale-dependent  time  interval  (conceptually  similar  to 
the  method  applied  in  [85]).  Consistent  with  the  notion  that  element  failure  corresponds  to  traversal 
of  a  macroscale  “train”  of  coalescing  microcracks  across  the  element,  a  small  element  fails  sooner 
than  a  large  element,  and  hence,  the  model  keys  the  “time-to-failure”  to  the  element  size.  For  sim¬ 
plicity,  elements  with  large  aspect  ratios  are  avoided  so  that  the  time-to-failure  has  a  length  scale 
of  approximately  where  V  is  the  volume  of  the  finite  element  that  is  softening.  Not  only  will 
larger  elements  be  initialized  with  lower  median  strength,  this  straightforward  (and  computation¬ 
ally  efficient)  method  for  scaling  the  time-to-failure  ensures  that  larger  elements  might  still  require 
greater  fracture  energy,  depending  on  how  the  material  parameters  are  set.  Implementation  details, 
including  formulas  for  degrading  strength  and  stiffness  with  damage,  are  provided  in  [6]. 

The  time-to-failure  scaling  is  insensitive  to  the  loading  rate.  Accordingly,  a  specimen  loaded 
slowly  will  seem  to  lose  strength  suddenly  because  the  time-to-failure  will  be  very  small  relative  to 
the  total  loading  time.  If  the  same  specimen  is  loaded  very  rapidly  to  the  same  strain  level,  it  will 
lose  strength  over  a  larger  fraction  of  the  total  loading  interval  (because  a  fixed  amount  of  time  must 
pass  before  the  stiffness  and  strength  degradations  fully  take  effect).  Thus,  when  loaded  rapidly, 
our  model  will  appear  to  have  a  finite  plastic  strain-to-failure  (which  is  similar  in  character  to  the 
Johnson-Holmquist  class  of  ceramic  models  [45, 48-50]  whenever  those  models  are  compared  with 
ours  al  a  given  loading  rate,  but  a  key  distinction  is  that  our  model’s  use  of  time-to-failure  provides 
different  apparent  strain-to-failure  at  different  loading  rates). 


Figure  18.  Laboratory  observations  of  the  effect  of  impact  speed  on  the  number  of  radial  cracks  [70].  Shown 
are  B4C  samples  impacted  at  speeds  ranging  from  100  to  400  m/s  from  left  to  right. 


Figure  19.  Realistic  trends  in  damage  patterns  when  including  spatial  variability,  scale  effects,  and  time-to- 
failure.  Changing  the  WC  impactor  speed  (from  200  to  500  m/s  from  left  to  right,  using  the  same  mesh  and 
property  seeding)  produces  a  corresponding  increased  number  of  radial  cracks  in  simulations  of  SiC  targets. 
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Because  our  model  requires  a  fixed  amount  of  time  to  pass  before  loss  of  strength  occurs,  the 
higher  strain  rate  caused  by  higher  impact  speed  allows  more  material  to  reach  the  failure  threshold 
before  release  waves  (or  Mott  waves — see  also  the  related  “event  horizon”  [76])  can  arrive  from  the 
radial  cracks  that  formed  slightly  sooner  at  weaker  points.  The  well-known  experimental  observation 
that  the  number  of  radial  cracks  in  axisymmetric  loading  increases  with  loading  rate  [86,  87]  (see, 
for  example,  observations  for  B4C  in  Figure  18)  is  therefore  easily  predicted  with  this  time-to-failure 
softening  algorithm  because  it  essentially  allows  a  larger  number  of  flaws  to  become  activated  at 
high  rates  (see  similar  trends  in  Figure  19  for  simulations  of  SiC).  Consistent  with  [58],  this  higher 
probability  of  flaw  activation  at  higher  loading  rates  was  also  mentioned  by  Brannon  and  Gowen 
[31]  as  being  a  simple  case  of  a  larger  portion  of  the  3D  stress  state  (as  depicted  in  our  Figure  4) 
falling  transiently  farther  above  each  crack’s  size-dependent  failure  line,  hence  naturally  activating 
larger  populations  of  both  flaw  sizes  and  orientations. 

8.  NECESSITY  OE  BOTH  STRENGTH  VARIABILITY  AND  TIME-TO-FAILURE  SCALING 

As  seen  in  Eigure  20,  imposing  both  strength  variability  and  time-to-failure  scaling  significantly 
reduces  mesh  dependencies  of  damaged  zones,  crack  patterns,  and  rate  of  crack  propagation  for  the 
dynamic  indentation  of  ceramics.  The  images  in  Figure  20  depict  simulations  of  a  6.35  mm  diameter 
tungsten-carbide  sphere  launched  at  375  m/s  into  a  25.4  mm  diameter  SiC-N  ceramic  cylinder, 
described  in  [7].  Each  set  of  images  in  Figure  20(a)-(d)  contains  a  low  resolution  simulation  next  to 
successively  higher  resolutions.  All  simulation  parameters  were  kept  the  same,  except  for  the  mesh 
resolution  and  the  feature  being  explored  in  the  model  (which  was  changed  for  each  subset). 

Figure  20(a)  contains  results  for  a  conventional  model  (Kayenta)  calibrated  to  the  median  of  the 
experimental  data  and  ignoring  size  effects.  Here,  the  high-pressure  limit  surface  parameter  Y  was 
calibrated  to  the  data  used  in  [49],  the  tensile  strength  parameter  h  was  calibrated  to  buttonhead  ten¬ 
sion  results  in  [88],  and  the  shear  strength  parameter  s  was  calibrated  to  triaxial  compression  results 
in  [11].  To  create  an  objective  comparison  between  different  features  of  the  model,  the  time-to- 
failure  was  computed  for  the  average  element  size  of  the  mesh  in  the  medium  resolution  simulation 
for  the  fully  enhanced  model  in  Figure  20(d)  (the  middle  image)  and  applied  as  a  constant  time-to- 
failure  for  the  conventional  model.  In  Figure  20(a),  the  number  of  radial  cracks  and  especially  the 
depth  of  penetration  depend  on  the  element  size  in  the  mesh.  For  the  simulations  in  Figure  20(b),  the 
time-to-failure  was  scaled  according  to  element  size  as  discussed  in  Section  7,  while  other  parame¬ 
ters  were  the  same  as  in  the  conventional  model  in  Figure  20(a).  Note  the  reduction  in  the  number 


(a)  Conventional  model  (no  variability  or  scaling) 


(b)  Time-to-failure  scaling  without  variability 


(c)  Variability  without  time-to-failure  scaling 


(d)  Variability  and  time-to-failure  scaling 


Figure  20.  Influence  of  strength  variability  and  time-to-failure  scaling  in  dynamic  indentation.  Subfigures 
each  show  three  results  at  the  same  moment  in  time  on  progressively  finer  mesh  resolutions. 
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of  radial  cracks  in  the  high-resolution  simulation  on  the  right  due  to  the  lower  time-to-failure  for 
the  smaller  elements  in  that  mesh.  For  the  simulations  in  Figure  20(c),  a  Weibull  variability  of 
strength,  including  the  volume  scaling  term  in  Equation  (7),  was  applied  to  the  h  parameter  in 
the  Kayenta  model,  as  discussed  in  Section  4.  The  other  parameters  (including  constant  time-to- 
failure)  were  the  same  as  in  the  conventional  model  in  Figure  20(a).  Note  the  significant  effect  of 
including  spatial  variability.  The  size  scaling  effects  have  resulted  in  a  significantly  higher  median 
element  strength  (and  hence  more  realistic  ballistic  performance)  while  the  variability  has,  in  prin¬ 
ciple  (modulo  issues  with  the  progression  of  damage  discussed  in  Section  7  that  can  be  encountered 
in  simulations  of  calibration  experiments),  preserved  the  match  of  the  model  to  the  suite  of  calibra¬ 
tion  experiments  discussed  in  Section  5.  Also,  using  a  size-independent  time-to-failure  causes  the 
horizontal  lateral  cracking  (on  the  perimeter  toward  the  bottom  of  the  target)  to  be  located  at  suc¬ 
cessively  higher  elevations  as  the  mesh  resolution  increases  (contrast  this  to  the  location  of  these 
cracks  in  Figure  20(d)). 

The  model  in  Figure  20(d)  includes  both  the  spatial  variability  of  strength  discussed  in 
Figure  20(c)  and  the  time-to-failure  scaling  discussed  for  Figure  20(b).  The  two  effects  work 
together,  resulting  in  a  model  with  substantially  reduced  mesh  dependence,  both  in  the  amount  of 
damage  induced,  and  in  the  rate  of  growth  of  damage.  The  approximate  number  of  radial  cracks  and 
location  of  lateral  cracking  are  similar  between  the  meshes.  Mesh  dependence  has  been  reduced,  but 
not  eliminated.  Significant  work  is  still  required  for  truly  predictive  ballistic  simulations  of  ceram¬ 
ics,  but  spatial  variability,  size  scaling,  and  time-to-failure  scaling  represent  technical  advances  in 
comparison  with  predecessor  theories.  Moreover,  these  features  can  be  added  without  significant 
increases  in  computational  (or  implementation/development)  overhead. 

9.  CONCLUSIONS 


Enhancing  an  existing  conventional  deterministic  damage  model  to  support  spatial  variability 
requires  (in  principle)  no  revisions  of  the  existing  constitutive  model’s  source  code.  Code  changes 
are  applied  at  the  host-code  level,  where  fields  are  allocated  (during  initialization)  for  spatially 
varying  material  parameters.  Eor  example,  each  element  in  a  simulation  might  have  its  own  ran¬ 
domly  assigned  scale-dependent  strength.  These  perturbed  parameters  are  assigned  values  by 
applying  whatever  statistical  distribution  is  evident  in  laboratory  testing  of  finite-sized  specimens 
but  with  care  given  to  anticipate  parameter  correlations  required  to  produce  admissible  material 
realizations  (e.g.,  to  avoid  inadvertently  generating  non-convex  strength  limit  surfaces). 

If  a  material  parameter  rj  happens  to  be  observed  in  the  laboratory  to  be  Weibull  distributed  and  if 
it  also  has  corresponding  Weibull  scale  effects,  then  the  spatially  variable  and  size-dependent  value 
assigned  to  a  hnite  element  of  volume  V  is  found  by  computing  a  random  number  R  and  then  saving 


r)  —  r] 


VlnR 

Uln(l/2) 


1  /m 


(18) 


into  a  field  array,  which  then  provides  each  finite  element  with  its  own  perturbed  material  property 
rj.  Here,  rj  is  the  median  value  of  rj  observed  in  the  laboratory  using  a  specimen  of  volume  V,  and 
m  is  the  Weibull  modulus  quantifying  observed  variability  in  the  property.  Appropriate  revisions 
naturally  apply  if  the  parameter  is  not  Weibull  distributed. 

Given  the  overarching  goal  to  introduce  variability  into  an  already  available  deterministic  dam¬ 
age  model,  each  material  property  rj  that  is  already  available  in  existing  property  databases  is  merely 
re-interpreted  as  the  median  about  which  uncertainty  and  scale  effects  are  imposed.  Thus,  the  only 
new  inputs  are  V  and  (if  Weibull-distributed)  m,  both  of  which  can  be  determined  by  reinspect¬ 
ing  the  laboratory  data  originally  used  to  determine  rj.  “Data  delocalization”  schemes  [38,  74-76] 
are  needed  to  set  the  effective  specimen  volume  V  for  any  strength  experiments  involving  non- 
uniform  stress  fields,  whereas  “relocalization”  [32]  essentially  re-introduces  stress  variation  that  is 
not  resolved  on  low-order  finite  elements. 

If  the  value  of  a  property  rj  was  determined  from  only  a  single  experiment,  then  the  experiment 
would  need  to  be  repeated  enough  times  to  quantify  the  scale-dependent  central  tendency  (e.g.. 
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median)  and  the  aleatory  uncertainty  (e.g.,  standard  deviation).  At  least  20  (and  preferably  more 
than  50)  such  tests  are  recommended  for  each  loading  mode  (shear,  tension,  compression,  etc.),  but 
indirect  methods,  such  as  “line  VISAR”  in  spall  testing  seem  promising  to  quantify  variability  in  a 
single  experiment  [53,  89]. 

Overall,  parameterizing  material  models  from  observed  aleatory  uncertainty  and  associated  scale 
effects  shows  promise  for  enhancing  the  predictive  capability  of  numerical  methods  of  fracture  and 
fragmentation  while  still  keeping  computational  and  development  costs  manageable.  An  ordinary 
deterministic  limit  surface  is  effectively  replaced  with  a  family  of  level  sets  (called  aleatory  quan¬ 
tile  surfaces)  associated  with  a  notional  failure  probability  function.  Various  scalar-based  strengths 
(e.g.,  buttonhead  tension,  triaxial  compression,  spall,  simple  shear,  Brazilian  indirect  tension,  etc.) 
each  have  different  distributions,  yet  they  are  all  simultaneously  modeled  by  reinterpreting  them 
as  directed  derivatives  of  the  single,  unihed  tensor-based  quantile  function.  The  tensorial  nature  of 
stress  is  therefore  directly  and  conveniently  accounted  for,  giving  non-Weibull,  pressure-dependent, 
and  generally  loading-mode-dependent  strength  distributions  based  on  laboratory  evidence.  A  scale- 
dependent  time-to-failure  algorithm  for  degradation  of  stiffness  (and  strength)  is  shown  to  not  only 
alleviate  mesh  sensitivity  in  dynamic  indentation  but  it  also  naturally  ensures  activation  of  larger 
populations  of  flaws  during  higher-rate  loading,  automatically  giving  more  fragments  with  increased 
loading  rate. 
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